A comparative analysis of the sensitivity and BOLD contamination of the VASO response at 3 Tesla: ME-DEPICTING vs. ME-EPI readouts^†

VAscular Space Occupancy (VASO) functional magnetic resonance imaging (fMRI) provides indirect measures of changes in cerebral blood volume (CBV; in units of ml blood per ml of tissue) by nulling all blood signal at the time of acquisition (Lu et al., 2003), such that the VASO signal is related to the CBV according to (Huber, Ivanov, et al., 2014):

Since its inception in 2003, the technique has evolved into many variants, including Multiple Acquisitions with Global Inversion Cycling (MAGIC) VASO, which allowed for multi-slice and whole brain imaging (Lu, Van Zijl, et al., 2004; Scouten & Constable, 2007); inflow-based VASO (iVASO) for arterial and arteriolar quantitative CBV measurements (Hua et al., 2011); or Slab-selective Inversion (SI) VASO (Jin & Kim, 2008) and Slice-Saturation Slab-Inversion (SS-SI) VASO (Huber, Ivanov, et al., 2014) for improved sensitivity. Of these, SS-SI-VASO has caught the most traction in the past few years due to its applicability in high-resolution high-field (≥7 T) layer-fMRI studies (Huber et al., 2015). The higher sensitivity and temporal resolution of SS-SI-VASO can also be exploited at lower field strengths. Corresponding implementations, however, were delayed due to the lower demand for layer-fMRI studies at 3 T. There are at present only three published SS-SI-VASO studies at 3 T (Guidi et al., 2023; Huber, Kronbichler, et al., 2023; Knudsen et al., 2023).

Alternatives to blood oxygenation level-dependent (BOLD) fMRI techniques, such as VASO fMRI, also serve as proxies for neuronal activation and help supplement a more holistic understanding of BOLD fMRI. The irony, however, is in the fact that at non-zero echo time (TE), these non-BOLD methods themselves suffer from contamination by the BOLD response (Hetzer et al., 2011; Huber et al., 2019). The functional sensitivity and accuracy of these measurements could then benefit greatly from short-TE readouts. One such readout is the Double-shot Echo Planar Imaging with Center-out Trajectories and Intrinsic NaviGation or DEPICTING (Hetzer et al., 2011). The multi-echo (ME) version of DEPICTING with very short first TE (⁠$TE1$⁠) and inter-echo time was recently found to substantially reduce BOLD contamination in pseudo-Continuous Arterial Spin Labeling (pCASL) (Alsop et al., 2015; Dai et al., 2008; Lorenz et al., 2018) measurements of cerebral blood flow (CBF) changes, while providing reliable measures of the simultaneous BOLD response (Devi et al., 2022). SS-SI-VASO generally takes care of the inherent BOLD contamination at non-zero TE by employing a dynamic division strategy, wherein the blood-nulled image is divided by its consecutive non-nulled (BOLD-weighted) control image. A complete correction for BOLD contamination, however, relies on a number of factors, most importantly being the assumption of a similar BOLD contribution in both the nulling and control condition, that is, extravascular (EV) BOLD contributions being equivalent to those of the combined extravascular plus intravascular (IV) BOLD response in the non-nulled parenchyma (Huber, Ivanov, et al., 2014).

In the present study, the feasibility of ME-DEPICTING as a potential readout for SS-SI-VASO at 3 T was investigated. Owing to its shorter TEs, a higher functional sensitivity is expected, as compared to the ME variant of the commonly used Echo Planar Imaging (EPI). All assessments of ME-DEPICTING-based SS-SI-VASO were, hence, contrasted to those of ME-EPI-based SS-SI-VASO. Apart from the provision for simultaneous BOLD and CBV measurements, the ME readouts also allowed for a comparison of the BOLD-correction strategies. In particular, a correction provided by ME extrapolation to $TE=0$ was tested against the correction provided by the original dynamic-division strategy.

Sixteen healthy volunteers (30 ± 5 years, 9 female) gave written informed consent before undergoing the experiments that had been approved by the Ethics Committee at the Medical Faculty of Leipzig University. All participants were right-handed and had normal or corrected-to-normal vision.

A full-field 8-Hz flickering black-and-white radial checkerboard was used for visual stimulation. Each functional cycle lasted 20 repetitions, starting with a rest block, which consisted of a blank gray screen of 12 repetitions (24 s), followed by the task block of 8 repetitions (16 s). The paradigm was programmed using Presentation (v17.2, Neurobehavioral Systems, Berkeley, CA, USA). A central, colored fixation point was present throughout the experiment. Subjects were instructed to focus on this dot and press a button whenever it changed color. Their attention was monitored by visually tracking their responses.

SS-SI-VASO was implemented on a 3-T MAGNETOM Skyra^fit scanner (Siemens Healthineers, Erlangen, Germany) equipped with a 32-channel receive head coil. Data were obtained from 10 slices of 4 mm thickness (no slice gap, nominal in-plane resolution 3 mm × 3 mm, field of view 192 mm, matrix 64 × 64, bandwidth 2232 Hz/Px) located along the calcarine sulcus and acquired in descending slice order. Selective inversion for the VASO scans (nulling condition) was achieved over a 34-cm slab centered at the middle of the slice package through a hyperbolic secant inversion radiofrequency (RF) pulse of 8 ms duration, bandwidth-time product of 10, and RF field peak amplitude of $B1$ = 14 μT. The blood-nulling condition, corresponding to an inversion time, TI = 1153 ms, for an assumed longitudinal relaxation time of blood of $T1,b=1664ms$ (Lu, Clingman, et al., 2004), was fulfilled for the 6th slice. A schematic of the SS-SI-VASO implementation is given in Figure 1 along with the sequence diagrams for ME-DEPICTING and ME-EPI. The specifications for the two two-dimensional (2D) ME readout modules were as follows:

• ME-DEPICTING: $TE1$/$TE2$/$TE3$ = 1.7 ms/10.7 ms/19.7 ms; repetition time, $TR$ = 2 s; GRAPPA (GeneRalized Autocalibrating Partial Parallel Acquisition) factor 2.

• ME-EPI: $TE1$/$TE2$/$TE3$ = 7.5 ms/20.7 ms/33.9 ms; $TR$ = 2 s; GRAPPA factor 2; partial-Fourier factor 6/8.

Functional runs consisted of ten functional cycles (200 repetitions, i.e., 100 nulling/control condition pairs; effective TR, $TReff$ = 4 s) and were of approximately 7 min duration. ME-EPI and ME-DEPICTING acquisitions were recorded during the same session. Their order was shuffled across participants to avoid primacy bias.

In addition to these two functional runs, two auxiliary functional runs of similar duration were recorded in one of the subjects, with ME-EPI and ME-DEPICTING specifications different from those stated above. They were:

• (i)

  ME-EPI with 5 echoes: $TE1$/$TE2$/$TE3$/$TE4$/$TE5$ = 10 ms/25.3 ms/40.6 ms/55.9 ms/71.2 ms; $TR$ = 2 s; GRAPPA factor 2; partial-Fourier factor 7/8.

• (ii)

  High-resolution ME-DEPICTING: 15 slices, 1.7 mm isotropic nominal resolution; $TE1$/$TE2$/$TE3$ = 2.8 ms/15.1 ms/27.4 ms; $TR$ = 2 s; GRAPPA factor 2.

For co-registration purposes, a 2D spoiled gradient-recalled echo (GRE) scan (1.5 mm nominal in-plane resolution; $TE$ = 3.1 ms; $TR$ = 1300 ms; flip angle 90°) was obtained at the start of each session, with slice geometry identical to the functional scans.

Data were preprocessed and analyzed using Statistical Parametric Mapping (SPM12; Wellcome Trust Centre for Neuroimaging, UCL, London, UK) implemented in MATLAB 2022b (The MathWorks, Natick, MA, USA), FMRIB Software Library (FSL, (Jenkinson et al., 2012)), and additional scripts written in Interactive Data Language (IDL 8.1, Exelis Visual Information Solutions, Boulder, CO, USA).

All functional data were preprocessed in an identical manner. The data were split into VASO and BOLD datasets using ‘fslsplit’ and ‘fslmerge’. The series acquired at $TE1$ of each dataset was realigned, and the resulting realignment parameters applied to the corresponding remaining echoes. To allow for a meaningful comparison of functional sensitivity between the two readouts, only participants who exhibited minimal motion, defined as translation <1 mm and rotation <0.015 radians in both scans, were included in the analysis. The mean correlations between the realignment parameters of the VASO and BOLD scans in these subjects were also evaluated. The realigned time series were temporally high-pass filtered at a cut-off frequency equal to two effective functional cycles [1/(20 $×TReff$₎] and 3D-Gaussian filtered at a full width at half maximum (FWHM) equal to the nominal voxel size (3 mm or 1.7 mm). Extrapolated ME signals at TE = 0, $Sbn(TE→0)$ (please refer¹) and $Sctr(TE→0)$⁠, as well as the effective transverse relaxation rates, $R2,bn*$ and $R2,ctr*$⁠, of the VASO and BOLD datasets, respectively, were then extracted from the preprocessed ME image volumes acquired at $TEi$ with voxel intensity $S(TEi)$ via linear regression (LINFIT) of the expression $lnS(TEi)=−TEi×R2*+lnS(TE→0)$⁠.

Percent VASO signal changes were then evaluated from: (i) uncorrected first echo of the blood-nulled dataset, $Sbn(TE1)$⁠, (ii) blood-nulled data of all echoes corrected by dynamic division, such that $Sdd(TEi)$$=Sbn(TEi)​/​Sctr(TEi)$⁠, and (iii) the ME-extrapolated intercept, $Sbn(TE→0)$⁠. The simultaneously acquired BOLD signal was evaluated (i) in terms of $ΔR2*$ from the $R2,bn*$ and $R2,ctr*$ regression results, and (ii) in terms of percent BOLD signal changes from the weighted summation of all echoes, $Sctr(sum)$ (Poser et al., 2006; Posse et al., 1999). The weights were computed from $Sctr(TEi)$ and the underlying BOLD model according to the fitted $R2,ctr*$ of each voxel (see Eq. 6 of Poser et al., 2006).

A general linear model (GLM) with the canonical hemodynamic response function (double gamma function) was implemented in IDL for statistical analyses and applied to all data.

Activation-related percent signal changes were quantified based on the resulting $β$ coefficients, such that $δs=[(​β1−β0)​/​β0]100%$⁠. Percent VASO signal changes were converted into relative changes in cerebral blood volume by assuming a resting value of $CBVrest$ = 0.05 ml/ml and negligible contribution from cerebrospinal fluid (CSF) (Scouten & Constable, 2007). The imperfect blood nulling in the remaining slices was then compensated based on Eq. 2b therein, rewritten here as:

with $Q=(1−Mz,IV(TI)​/​Mz,EV(TI))−1$ for non-nulled slices; $Q=1$ for the nulled slice; and water densities of blood, $ρIV=$ 0.86 ml water/ml blood, and parenchyma, $ρEV+IV=$ 0.89 ml water/ml parenchyma. $Mz,IV(TI)$ and $Mz,EV(TI)$ are the longitudinal magnetizations per unit of water in blood and parenchyma, respectively, at the time of acquisition; and $ΔCBV$ is the absolute $CBV$ change in ml blood/ml of parenchyma. For calculation of the quantitative rate changes $ΔR2*$ in units of s⁻¹, $R2*,rest$ was taken directly from the estimated $β0$ parameter of the respective $R2*$ time series.

Significant regions of VASO and BOLD activation were identified based on a voxel-based significance threshold of $p<10−6$⁠, except for the high-resolution data, which owing to the lower signal to noise ratio (SNR) was thresholded at $p<10−4$⁠. Slices cropped during the realignment process due to inter-session motion were disregarded from all scans of the particular session.

The sensitivity of the VASO response was evaluated from $Sbn(TE1)$⁠, $Sdd(TEi)$ and $Sbn(TE→0)$⁠, while that of the BOLD response was evaluated from $R2,ctr*$ and $Sctr(sum)$⁠. The comparison of $Sdd(TEi)$ and $Sbn(TE→0)$ VASO responses, inadvertently, allowed a comparison of the two BOLD correction strategies. The functional sensitivities were assessed in terms of (i) number of suprathreshold voxels, and (ii) the more reliable metric of temporal CNR (Geissler et al., 2007). The latter was evaluated from the common region of interest (ROI) between the ‘VASO activation’ obtained from the ME-extrapolated $Sbn(TE→0)$ and the ‘BOLD activation’ taken from the $R2,ctr*$ images of both readouts. Comparisons of relative VASO or BOLD signal changes across TEs were also based on this ROI.

The parenchymal and extravascular BOLD contributions were estimated to evaluate the BOLD correction by dynamic division, which relies on the assumption $T2,EV+IV*$$≈$$T2,EV*$⁠. The intravascular BOLD contribution was then estimated in terms of $ΔR2*$⁠, written as $ΔR2,IV*$⁠. These analyses were based on two ROIs in the nulled slice (i.e., the 6th/8th slice for standard-/high-resolution data):

• (i)

  ROI-1 was based on common BOLD activation in the $R2,ctr*$ data of the two readouts. $ΔR2*$ evaluated from $R2,ctr*$ was expected to represent parenchymal $ΔR2*$ and is henceforth written as Δ $R2,EV+IV*$⁠.

• (ii)

  ROI-2 was based on common BOLD activation in the $R2,bn*$ data of the two readouts. Δ $R2*$ evaluated from $R2,bn*$ was expected to represent extravascular $ΔR2*$ and is henceforth written as Δ $R2,EV*$⁠.

The related effective transverse relaxation times during the resting condition, $T2,EV+IV*,rest$ and $T2,EV*,rest$⁠, were also extracted. Prior to the extraction of the blood-nulled slice from the visual cortex masks (-fslslice), the multi-slice mask was multiplied by a brain mask that had been extracted from the structural image (-bet with -f 0.2) (Smith, 2002) and linearly registered with the $Sbn(TE→0)$ images of each readout (-flirt with -applyxfm, -usesqform & -noresampblur) (Jenkinson et al., 2002; Jenkinson & Smith, 2001). The relative extravascular rate change was evaluated from ROI-1 and ROI-2, such that

The effective intravascular BOLD contribution of both readouts was estimated solely from ROI-2 according to:

An expression for $E2,IVact≡exp(−TE×R2,IV*,act)$ is derived in the Appendix A1 (Eq. A6). Experimental values of $T2,EV*,rest$ and $T2,EV*,act$ were taken from mean values in ROI-2. Similarly, $CBVact$ was calculated from $δsVASO,bn(TE→0)$ and the VASO signal from the dynamically divided middle echo $δsVASO,dd(TEmid)$ with $TEmid= TE2$ = 10.7 ms and 20.7 ms for ME-DEPICTING and ME-EPI, respectively, and $TEmid=TE3$ = 40.6 ms for the ME-EPI with 5 echoes. The other relevant parameters had the following assumed values: $T1,EV$ = 1330 ms (Wansapura et al., 1999); relative proton densities, $ρIV$ = 0.87 ml water/ml blood; $ρEV+IV$ = 0.89 ml water/ml parenchyma (Donahue et al., 2009; Herscovitch & Raichle, 1985; Lu et al., 2002). $ΔR2,IV*$ was estimated for $T2,IV*,rest$ varying from 15 ms to 40 ms in steps of 5 ms, that is from a more venous to a more arterial regime (Zhao et al., 2007). The influence of resolution and longer TE on $ΔR2,IV*$ and its dependence on $T2,IV*,rest$ was also examined for a single subject with standard versus high-resolution ME-DEPICTING data and 3-echoes versus 5-echoes ME-EPI data.

Four participants were excluded from the evaluation based on the motion criterion and another one due to minimal activation. Supplementary Figure S1 shows examples of plots of the translational and rotational displacements of an excluded and an included subject. The mean correlation between the realignment parameters of VASO and control scans over the remaining 11 participants averaged at 0.86 ± 0.07 and 0.86 ± 0.09 for EPI and DEPICTING scans, respectively (Supplementary Table S1).

Significant VASO and BOLD activations were identified in all remaining participants. Examples are presented in Figure 2. The higher sensitivity of $Sbn(TE1)$ of DEPICTING at the subject level (Fig. 2A) was confirmed for the group-averaged data (Fig. 3A). The number of significant voxels showing VASO, BOLD and their common activation (VASO $∩$ BOLD, from $Sbn(TE→0)$ and $R2,ctr*,$ respectively) for each participant can be assessed from Supplementary Table S2. The temporal VASO CNR values are shown in Figure 3B and Supplementary Table S3. Interestingly, after extrapolation to zero TE, the activated areas of both readouts became rather similar and the CNR was almost equal.

It is to be noted that even though the ROI selection relied on BOLD activation from $R2,ctr*$ (Supplementary Table S2), the sizes of the ROI were similar even when the significant BOLD area was based on $Sctr(sum)$ (not shown). The sensitivity difference for the detection of BOLD signal changes was found to be negligible (~7%) between the $Sctr(sum)$ and $R2,ctr*$ for both readouts (Fig. 3C and Supplementary Table S4). The CNR of BOLD data obtained from each $Sctr(TEi)$ for the two readouts from the ROI is plotted in Figure 3D.

Investigation into the startling gain in VASO sensitivity of the dynamic division-corrected data revealed further increases in both the activation area (Supplementary Table S2) and VASO CNR (Supplementary Table S3) for $Sdd(TE2$⁠) and $Sdd(TE3$⁠) for both readouts. Significantly activated areas in two arbitrary participants are presented in Figure 4 along with corresponding maps obtained from $Sbn(TE→0)$ and $Sbn(TE1)$⁠. However, in case of the dynamic-division method, along with the increase in activation area with TE, the intensity of the VASO signal change, $δsVASO,dd$ was also found to increase with TE. The increase of $δsVASO,dd$ with $TE$ is also evident in their cycle-averaged time courses (green line with filled squares in Fig. 5).

The time courses in Figure 5 were extracted separately for each readout from a common ROI of VASO activation in the nulled slices of ME-DEPICTING and ME-EPI data. As expected, $δsVASO,bn(TEi)$ switches from a negative to a positive response with increasing $TE$ for both readouts, while $δsBOLD,ctr(TEi)$ grows with $TE$⁠, as evidenced by the amplitudes of the primary signals and post-stimulus undershoots. Interestingly,$δsVASO,dd(TEi)$ can be seen mirroring the corresponding BOLD signals. This applies to the evolution of both the signal intensities and the shapes of the time courses as a function of $TE$ for both readouts. Another finding is that $δsVASO,dd(TEi)$ agreed rather well with the difference between the uncorrected VASO response during the nulling condition and the BOLD response during the control condition, $δsVASO,bn(TEi)−δsBOLD,ctr(TEi)$⁠.

Figures 6 and 7 further demonstrate the TE dependence of the dynamic division-corrected VASO data for both readouts. With a near-perfect linear increase (⁠$r2≥0.97$⁠), the relation of $|δsVASO,dd|$ bears a striking resemblance to the TE dependence of $δsBOLD,ctr$⁠. The slopes of the linear fits of $|δsVASO,dd|$ and $δsBOLD,ctr$ averaged across subjects were found to be comparable between the two readouts.

Furthermore, the $TE$ dependence of $δsVASO,dd$ was also identified for the auxiliary single-subject data with 5 echoes (Fig. 8). A reduction in the extent of activation can be seen with later echoes (Fig. 8A). The amplitude of the primary signal of VASO responses grows rapidly at smaller TEs before appearing to stabilize at $TE$ ≥ 25.3 ms. A shift in the shape of the time courses can be seen for $δsVASO,dd(TEi)$⁠; specifically, in a faster return to baseline compared to $δsVASO,bn(TE1)$ and $δsVASO,bn(TE→0)$⁠. An increase in the amplitude of the post-stimulus overshoot is also evident (Fig. 8B). Figure 9 demonstrates a similar behavior of the dynamic division-corrected data (increase of $δsVASO,dd$ with $TE$⁠) for higher-resolution data. The shape of the signals from $δsVASO,dd(TE1)$ deviates again from that of uncorrected $δsVASO,bn(TE1)$ and ME-extrapolated $δsVASO,bn(TE→0)$⁠. The larger amplitude of VASO signal changes at higher resolution is also evident.

Averaged functional $ΔR2,EV*$ and $ΔR2,EV+IV*$ values resulted in estimates of the extravascular $ΔR2*$ fraction, over the parenchymal ROI-1 (146 ± 59 voxels), of $fEV$ = 55±15% and 50 ± 12% for ME-EPI and ME-DEPICTING, respectively. A substantial increase in this fraction was obtained in the smaller extravascular ROI-2 (45 ± 31 voxels): $fEV$ = 62 ± 14% and 60 ± 11%, respectively. The sizes of ROI-1 and ROI-2 for each participant and the corresponding $ΔR2,EV+IV*$⁠, $ΔR2,EV*$ and $fEV$ are provided in Table 1 for both readouts along with the $T2,EV+IV*,rest$ and $T2,EV*,rest$ values. $T2,EV*,rest$ was found to be slightly shorter than their $T2,EV+IV*,rest$ counterparts with both EPI and DEPICTING in both ROIs. An example of the voxel-wise distribution of $T2,EV+IV*,rest$⁠, $T2,EV*,rest$⁠, $ΔR2,EV+IV*$⁠, $ΔR2,EV*$ and the extravascular $ΔR2*$ fraction over ROI-1 (participant P16: average $fEV$ of 51% and 46% for EPI and DEPICTING, respectively) is shown in Figure 10.

Increased estimates of the extravascular $ΔR2*$ fraction were obtained from the single-subject high-resolution DEPICTING data with $fEV$ = 80% (ROI-1, 458 voxels) and 95% (ROI-2, 237 voxels). A somewhat higher extravascular $ΔR2*$ fraction was also obtained with the 5-echo EPI data. These results along with those reported in recent literature for the visual cortex are summarized in Table 2.

Table 3 lists the experimentally obtained data employed in the estimation of intravascular $ΔR2,IV*$ with Eq. 4. Subject-averaged values have been provided for the multi-subject main study with 3-echoes and standard resolution. Values, corresponding to the participant (P1) from this study have also been provided along with their values for the single-subject study with 5 echoes (ME-EPI) and higher resolution (ME-DEPICTING). All values were based on ROI-2. The $CBVact$ (⁠$CBVrest+ΔCBV$⁠) values were obtained by assuming $CBVrest$ = 0.05 ml/ml and are based on the $Sbn(TE→0)$ data. Subject-averaged $ΔCBV$ values of 0.0053 ml/ml and 0.0067 ml/ml corresponded to percent CBV changes of $δcbv$ = 13 ± 3 and 11 ± 4% for EPI and DEPICTING, respectively. The difference between the results obtained with both readouts was insignificant (paired two-tailed t-test, $p$ = 0.08). Similarly, $δcbv$ = 11% was obtained for the 5-echo ME-EPI while the high-resolution DEPICTING data yielded the highest $CBV$ change, $ΔCBV=$ 0.0135 ml/ml or $δcbv$ = 27%.

The results for the intravascular $ΔR2,IV*$ estimations are plotted in Figure 11 for assumed $T2,IV*,rest$ values ranging from 15 to 50 ms. A nonlinear relationship between $ΔR2,IV*$ and $T2,IV*,rest$ is evidenced in all the plots. A steeper decline of $|ΔR2,IV*|$ is observed in the region of venous $T2,IV*,rest$⁠. The differing $ΔR2,IV*$ contributions between the two readouts can also be seen converging at higher resting $T2*$ values (Fig. 11A). The within-sequence comparisons demonstrate a lower intravascular contribution (i.e., smaller $|ΔR2,IV*|$⁠) with higher spatial resolution (Fig. 11B) and longer $TE$ (Fig. 11C).

Readouts with short TE come with the promise of better sensitivity for non-BOLD contrasts (Hetzer et al., 2011; Huber et al., 2019) as BOLD contamination at $TE$ <2 ms is expected to be minimal. Our recent pCASL study at 3 T had proven this to be the case for functional CBF changes in the visual cortex (Devi et al., 2022). The area of activation obtained with pCASL-prepared ME-DEPICTING at $TE1$ = 1.7 ms exceeded that obtained with ME-EPI at $TE1$ = 8 ms by 40%. The present study confirms the same for the VASO contrast. The sensitivity for detecting uncorrected VASO signal changes at $TE1$ of the blood-nulled dataset, $Sbn(TE1=1.7 ms)$ was found to be substantially higher than that of ME-EPI, $Sbn(TE1=7.5 ms)$⁠. However, comparing the functional sensitivity of data BOLD-corrected by ME extrapolation to $TE=0$ (Fig. 3B), an almost equivalent VASO CNR was obtained for both readouts. This differs from the results of the pCASL study, wherein a higher CBF-CNR was obtained for the ME-DEPICTING acquisitions extrapolated to zero $TE$ [Fig. 6C and 6E in Devi et al. (2022)] with approximately the same TEs as employed in the current work. Also contrary to the VASO results was a higher CNR for CBF changes obtained with $TE1$ data for both readouts compared to the corresponding $S0$ data. An explanation for the improvement of the quality of the fitted data at $TE→0$ in SS-SI-VASO could be the fact that the data were first split into blood-nulling and control images prior to extrapolation, whereas the pCASL data were fit as a whole, with the fluctuations between the control and label data, possibly degrading the quality of the resulting $S0$ estimate. Nonetheless, as with the pCASL study, the sensitivity for the BOLD response obtained from the $R2,ctr*$ and $Sctr(Sum)$ data was found to be comparable for both readouts (Fig. 3C and Supplementary Table S4). For instance, $R2,ctr*$ yields 1500 ± 448 suprathreshold voxels with ME-DEPICTING vs. 1420 ± 439 voxels with ME-EPI (⁠$p$ = 0.60; CNR 3.2 $±$ 0.9 vs. 3.7 $±$ 1.1, $p$ = 0.05). This finding, along with the prediction of a shorter $TE1$ providing a better approximation of CBV (Genois et al., 2021), makes ME-DEPICTING a promising substitute of the traditional EPI readout for the simultaneous measurement of CBV changes accompanying the BOLD response.

A much higher gain in VASO sensitivity was deduced for the $TE1$ data corrected for BOLD contributions by dynamic division compared to the correction by extrapolation to zero TE for both readouts. Further investigation into this remarkable improvement revealed a linear TE dependence of the BOLD-corrected $δsVASO,dd$ data (Figs. 6A, 6D, 7A and 7C). This TE dependence bears a close resemblance to that of the BOLD data, the difference being that the latter’s intercept is close to zero (Figs. 6C, 6F, 7B and 7D). This suggests a failure of the dynamic division strategy and, hence, the presence of residual BOLD contamination in the $δsVASO,dd$ data, as also hinted by the similar but mirrored shape of their time courses. The VASO measures obtained from extrapolated and even uncorrected data followed temporal dynamics expected for a functional CBV response: a slower return to baseline and less prominent post-stimulus transients (Lu et al., 2003; Mandeville et al., 1998). Dynamic division-corrected data, on the other hand, showed an earlier return to baseline and post-stimulus overshoots that increased in peak amplitude with TE (Figs. 5 and 8B). This finding mimics that of the BOLD response in ME acquisitions (Havlicek et al., 2017) and was even more evident in the ME-EPI data with extended echo train length (⁠$TE5$ = 71.2 ms) (Fig. 8B).

Relating VASO signal changes to the CBV response in SS-SI-VASO (Eq. 1) relies on the assumption of equivalent BOLD contributions in the blood-nulling and control condition (Huber, Ivanov, et al., 2014). This assumption holds in situations where the BOLD signal is mostly of extravascular origin. The extent of intravascular BOLD contributions to the BOLD signal, however, relies on a number of factors, including field strength (⁠$B0$⁠), TE, readout sequence, and diffusion weighting (Duong et al., 2003; Silvennoinen et al., 2003; Uludaǧ et al., 2009). The intravascular $ΔR2*$ fraction is expected to decrease with $B0$ with predicted intravascular BOLD contributions from the microvasculature of approximately 57%, 36%, 11%, and 5% at the respective field strengths of 1.5 T, 3 T, 4 T, and 4.7 T and $TE$ = $T2,EV*$ (Uludaǧ et al., 2009). Experimentally obtained extravascular $ΔR2*$ fractions are approximately 70% at 3 T (Donahue et al., 2011; Lu et al., 2003) and approximately 90% at 7 T (Cheng et al., 2015; Donahue et al., 2011) in the human visual cortex. The results from the current study are well within the range of what has been reported (Table 2). Our estimated $ΔR2,b*$ values (Fig. 11) agree also with previous in-vitro measurements in bovine blood of varying hematocrit (Hct) levels (⁠$ΔR2,b*$ = –8.16 s^–1,—14.3 s^–1, and—16.6 s^–1 for Hct = 0.21, 0.44, and 0.57, respectively), with the assumption of changes in blood oxygen saturation fraction, Y, from 0.61 to 0.73 (Zhao et al., 2007).

The observed linear dependence on $TE$ of the VASO data corrected by dynamic division align with simulations based on a vascular anatomical model (Genois et al., 2021), wherein the intravascular BOLD contributions were taken into account. The linear dependence of $δsVASO,dd$ on $TE$ ranging from 0 to 30 ms as depicted in Figure 5A of Genois et al. (2021) for 3-T data with an extravascular fraction of 72% resemble our plots in Figure 6A and 6D.

Our results, however, differ from those reported in the original SS-SI-VASO paper (Huber, Ivanov, et al., 2014). The 3-T results therein, with data from two participants, found $δsVASO,dd$ to be independent of TE (14 ms and 30 ms) in the visual cortex. This is remarkable as at a nominal resolution of 3 × 3 × 4 mm³, partial voluming effects are expected to be at play and introduce significant intravascular BOLD contributions. For instance, with the assumption of a 5% increase in oxygen metabolism, 57.5% increase in CBF, TE = 50 ms, and TR = 1s, Zhao et al. (2007) estimated an intravascular BOLD contribution of 13% from pure parenchyma, which increased to a substantial 42% for parenchyma contaminated by 2% veins. Consistently, a reduction in intravascular BOLD contributions from larger veins at higher resolution is evident in our single-subject data (Fig. 11B). Interestingly, despite an estimated 80–90% extravascular $ΔR2*$ fraction for this data, the intravascular fraction was sufficient to bring about a TE dependence of $δsVASO,dd$⁠.

A TE independence of 3-T ME-EPI data (12 ms and 48 ms) corrected by dynamic division has been reported in a recent layer-fMRI VASO study at sub-millimeter resolution (Huber, Kronbichler, et al., 2023), where contamination by larger veins is expected to be largely reduced. It should be noted that even in the absence of intravascular BOLD signal, any CBV change leads to a TE-dependent SS-SI-VASO signal induced by the TE dependence of the resting signals of corresponding (arterial or venous) blood compartments relative to the tissue signal. Consequently, an SS-SI-VASO signal increasing with $TE$ was also found at 7 T in a study employing a resolution of 1.5 mm (Fig. 5A in Huber, Goense, et al., 2014). Quantification of CBV changes during activation at $TE>0$ will always be governed by such biases. Under certain circumstances, however, these biases might cancel out for $TE>0$⁠, because such arterial and venous contributions tend to increase and decrease the SS-SI-VASO signal, respectively.

The echo-wise evaluation indicated a remarkable similarity of amplitude and transients of $δsVASO,dd$ and the difference signal, defined as $δsVASO,bn−δsBOLD,ctr$ (Fig. 5). Differences in the extravascular $ΔR2*$ fraction and, subsequently, in intravascular BOLD contaminations do not explain this finding. As shown in Appendix A2, the similarity of both quantities is easily explained by the relation $δsVASO,dd=(Sctrrest​/​Sctract)(δsVASO,bn−δsBOLD,ctr)$ (see Eq. A9). Because $Sctrrest≤Sctract$ holds (positive BOLD response) and given the small stimulus-related signal change (~1–2%), $|δsVASO,dd|$ is expected to be slightly smaller than $|δsVASO,bn−δsBOLD,ctr|$⁠. This relationship, however, could explain the stabilization of $δsVASO,dd$ at longer $TE$ in our 5-echo example. Since $δsBOLD,ctr$ contains an intravascular signal whereas $δsVASO,bn$ does not, their difference is expected to vary faster at shorter $TE$⁠, because at 3 T, the intravascular signal increases significantly more than the extravascular (Triantafyllou et al., 2011).

Differences between the EPI readouts employed in previous studies and the DEPICTING readout, such as different $T2*$-sensitivities of the specific k-space trajectories (Devi et al., 2022; Patzig et al., 2021), may contribute to the failure of the dynamic-division strategy in the present work. Although this aspect is likely to be of only secondary order because a TE dependence was introduced in our $δsVASO,dd$ data and corresponding $δcbf$ acquired with EPI as well (Fig. 6A, 6B, 6D and 6E), it cannot be ignored. As shown in Appendix A3, additional $T2*$-weighting due to the k-space trajectory would result in non-zero $δsctr(TE→0)$ and may be captured by an offset term $TEt$ leading to an effective echo time, $TEeff=TE+TEt$⁠. With the intercepts obtained in Figure 7B and 7D and the results from Table 2, Eq. A11 yields estimates of $TEt$ ≈ 0.96 ms and 0.35 ms for DEPICTING and EPI, respectively. The larger offset $TEt$ for DEPICTING is due to the different $T2*$-sensitivity of its k-space trajectory (Devi et al., 2022; Patzig et al., 2021). A significant deviation from the nominal TE results only at $TE1$ of DEPICTING (⁠$TEt$ ≈ 0.56 $TE1$⁠), but no relevant bias at longer TEs and for EPI. This effect leads to the slightly lower $δcbv$ value estimated from $Sbn(TE→0)$ of DEPICTING (11 ± 4%) compared to that of EPI (13 ± 3%) and also explains somewhat larger BOLD-like fluctuations in its extrapolated $δsctr$ signal (cf., participants P11 and P12 in Supplementary Fig. S2). The different zero-crossing TE of the $δsbn$ of the two readouts (Supplementary Fig. S3: approximately 11 ms and 17 ms for DEPICTING and EPI, respectively) could also be attributed to this difference in BOLD weighting at equivalent TEs.

A distinction lies in the TEs employed by earlier studies, which were comparatively longer and, hence, less likely to be significantly influenced by intravascular BOLD contributions. A consistent tendency could be the near-plateauing of $δsVASO,dd$ values at longer TE of the ME-EPI data, while the ME-DEPICTING data with shorter TEs appear to follow an increasing trend over the entire range (Fig. 6A and 6D). This is also evident in time courses of later echoes extracted from the nulled slice of ME-EPI in Figures 5 and 9B. Further ME studies might provide additional insight into this matter. Such studies are, however, presently rare, even with SS-SI-VASO being the current workhorse of CBV-based layer-fMRI studies. Apart from the work mentioned above, all other studies used single-echo acquisitions at 7 T with 15 ms ≤ $TE$ ≤ 28 ms (Beckett et al., 2020; de Oliveira et al., 2023; Faes et al., 2023; Guidi et al., 2016, 2020; Huber et al., 2015, 2021; Huber, Kassavetis, et al., 2023; Liu et al., 2022; Oliveira et al., 2021, 2022; Pizzuti et al., 2023) and with 27 ms at 3 T (Knudsen et al., 2023). We note that we do not expect relevant differences between the 2D readouts used here and those of the 2D Simultaneous Multi-Slice (SMS) EPI and 3D EPI variants (Huber et al., 2018), given that evaluations were exclusively based on the nulled slice.

Short TEs can also be achieved with spiral acquisitions. With recent developments in their implementation, artifact correction, and reconstruction, these can be expected to be a promising alternative to EPI sequences with Cartesian trajectories for non-BOLD fMRI (Glover, 2012; Kasper et al., 2022). A preliminary report of a two-fold improvement in temporal SNR (tSNR) compared to 3D-EPI was recently presented for VASO fMRI with spiral readouts at 7 T (Monreal-Madrigal et al., 2023).

SS-SI-VASO was implemented in the present study for TI/TR = 1153 ms/2000 ms. This fulfills the requirements of the original paper [see Fig. 1a of Huber, Ivanov, et al. (2014)] by choosing a TI shorter than the arterial transit time of blood and a ‘period III’ (i.e., $2TR+TI$⁠) that prevents contamination with blood that was inverted more than once. Our timing, however, does not null the CSF contribution to the signal. A higher extravascular signal intensity was preferred over CSF nulling by our choice of a shorter TR compared to that required for simultaneous nulling of blood and CSF contributions (TR = 2.75 s). Pilot experiments had revealed a more focal activation with CSF nulling. However, $δsVASO,bn(TE→0)$ values within a common ROI were found to be very similar. Nonetheless, partial volume effects due to CSF contribution cannot be disregarded in our results. The TI of 1153 ms was based on a blood $T1$ value of 1664 ms measured from bovine blood at 3 T (Lu, Clingman, et al., 2004), which is lower than recently published $T1$ values (~1800 ms) of human blood (Li et al., 2016, 2017). With our approximate slice acquisition time of 71 ms, our choice of a TI might has nulled the 7th rather than the 6th slice. This is expected to result in only minimal blood contributions and, hence in significant over- or underestimations of $CBV$ changes.

The ROI definition in the current study could also have been attributed to deviations from the original report. The voxel-based thresholding employed in the present work could bias the selection to voxels that showed the maximum percent signal changes and, hence, towards voxels with more partial voluming with blood vessels. On the other hand, the more widely employed cluster-based thresholding is sensitive to weaker and more diffuse signal changes, while suffering from lower spatial specificity (Woo et al., 2014).

The very short TE achieved with ME-DEPICTING was of primary interest for the current work and as such, 2D readouts were compared. While the error due to incompletely nulled slices was corrected for evaluating $ΔCBV$⁠, such corrections are not easily obtained for $ΔR2*$ values. Consequently, estimations of $ΔR2,EV*$⁠, and hence also $ΔR2,IV*$⁠, were limited to the single blood-nulled slice. Future work with 2D SMS-EPI (Huber et al., 2016) and SMS-DEPICTING or corresponding 3D readouts could solve this issue. At 7 T, 2D SMS-EPI was found to outperform 3D EPI (Poser et al., 2010) at lower resolutions due to lesser contributions from physiological noise (Huber et al., 2018).

Only a single pilot acquisition at higher resolution, which is likely to be less biased by intravascular BOLD contributions, and another one with more echoes have been provided here. Increasing the number of participants in each could further solidify these findings. Assessing the TE dependence may benefit from a further increase of the number of echoes. Extending the evaluation beyond the visual system, such as the motor cortex, is also warranted. The impact of intravascular BOLD contributions might also be investigated by combining 3-T experiments with those at 7 T, where the extravascular BOLD effect is expected to be dominant. Further experimentation is, indeed, needed to comment on suitable ranges of $TE$ and spatial resolutions, within which the dynamic division method would be feasible at the respective field strengths. This would help identify the impact of the dynamic-division correction in SS-SI-VASO-based layer fMRI studies.

Presently, in the absence of a ground truth, $δsVASO,bn(TE→0)$ and, subsequently, $δcbv$ values obtained by extrapolation to zero TE appear to deliver the most reliable estimations, especially, when quantification is the goal. However, the assumed mono-exponential $T2*$ decay with TE is likely to be too simplistic and may introduce some bias.

SS-SI-VASO was implemented at 3 T with two 2D readouts. The feasibility of the ME-DEPICTING sequence as a potentially advantageous implementation compared to ME-EPI was investigated, motivated by the shorter $TE1$ and inter-echo time capabilities of this double-shot EPI readout. Unwanted BOLD contamination was, in fact, drastically reduced at the very short $TE1$ of 1.7 ms. The extrapolation to zero TE, however, resulted in equivalent VASO sensitivity for the two readouts. Interestingly, the dynamic division approach was found to bring about a higher gain in VASO sensitivity for ME-EPI than for ME-DEPICTING. VASO signal changes corrected for BOLD with this strategy were, however, found to exhibit a distinct $TE$- dependence. This is probably due to the influence of intravascular BOLD contributions at 3 T. The correction for BOLD contamination by extrapolation to zero TE using multiple echoes does not rely on a negligible intravascular BOLD signal and is, therefore, recommended, if an overestimation of functional CBV changes is to be avoided. Consequently, in the absence of ME data, the short $TE$ of the DEPICTING readout still provides an alternative to conventional EPI if better sensitivity and accuracy of the VASO data are desired.

Pre-processed data for all echoes, data derived from multi-echo fitting and relevant scripts for BOLD correction by dynamic division and estimation of $ΔR2,IV*$ are available at https://osf.io/cg7sp/.

R.D.: Conceptualization, Data Acquisition, Methodology, Software, Formal analysis, Investigation, Writing—original draft, Writing—review & editing, and Visualization. J.L.: Methodology, Software. T.M.: Conceptualization, Methodology, Software, Investigation, Writing—original draft, and Writing—review & editing. H.E.M.: Conceptualization, Investigation, Writing—original draft, Writing—review & editing, Supervision, Project administration, and Funding acquisition.

The DEPICTING sequence is registered under Hetzer S, Mildner T, and Moeller H. 2014. Magnetic resonance imaging with improved imaging contrast. US Patent 8,664,954 B2, filed March 31, 2009, and issued March 4, 2014; European Patent 2 414 861 B1, filed March 31, 2009, and issued January 28, 2015.

We thank Dr. Laurentius (Renzo) Huber for helpful discussions on the implementation of SS-SI-VASO at 3 T. We also greatly appreciate Anke Kummer, Sylvie Neubert, and Nicole Pampus for their assistance with the volunteers and ascertaining smooth operation despite COVID regulations/restrictions. And we hope to memorialize Torsten Schlumm’s contribution to the development of the ME-DEPICTING project with our dedication.

A preliminary account of this work has been presented in the Proceedings of the 31st Annual Meeting of the ISMRM, London, UK, 2022. This work has been supported by the Max Planck Society and by the International Max Planck Research School on Neuroscience of Communication: Function, Structure, and Plasticity.

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00333.