Comparing hand movement rate dependence of cerebral blood volume and BOLD responses at 7T

Functional magnetic resonance imaging (fMRI) based on the Blood Oxygenation Level Dependent (BOLD) contrast takes advantage of the coupling between neuronal activity and the hemodynamics to allow a non-invasive localization of the neuronal activity. In general, fMRI experiments assume a linear relationship between neuronal activation and the observed hemodynamics. However, the relationship between BOLD responses, neuronal activity, and behaviour are often nonlinear. In addition, the nonlinearity between BOLD responses and behaviour may be related to neuronal process rather than a neurovascular uncoupling. Further, part of the nonlinearity may be driven by vascular nonlinearity effects in particular from large vessel contributions. fMRI based on cerebral blood volume (CBV), promises a higher microvascular specificity, potentially without vascular nonlinearity effects and reduced contamination of the large draining vessels compared to BOLD. In this study, we aimed to investigate differences in BOLD and VASO-CBV signal changes during a hand movement task over a broad range of movement rates. We used a double readout 3D-EPI sequence at 7T to simultaneously measure VASO-CBV and BOLD responses in the sensorimotor cortex. The measured BOLD and VASO-CBV responses increased very similarly in a nonlinear fashion, plateauing for movement rates larger than 1 Hz. Our findings show a tight relationship between BOLD and VASO-CBV responses, indicating that the overall interplay of CBV and BOLD responses are similar for the assessed range of movement rates. These results suggest that the observed nonlinearity of neuronal origin is already present in VASO-CBV measurements, and consequently shows relatively unchanged BOLD responses.


a b s t r a c t
Functional magnetic resonance imaging (fMRI) based on the Blood Oxygenation Level Dependent (BOLD) contrast takes advantage of the coupling between neuronal activity and the hemodynamics to allow a non-invasive localisation of the neuronal activity. In general, fMRI experiments assume a linear relationship between neuronal activation and the observed hemodynamics. However, the relationship between BOLD responses, neuronal activity, and behaviour are often nonlinear. In addition, the nonlinearity between BOLD responses and behaviour may be related to neuronal process rather than a neurovascular uncoupling. Further, part of the nonlinearity may be driven by vascular nonlinearity effects in particular from large vessel contributions. fMRI based on cerebral blood volume (CBV), promises a higher microvascular specificity, potentially without vascular nonlinearity effects and reduced contamination of the large draining vessels compared to BOLD. In this study, we aimed to investigate differences in BOLD and VASO-CBV signal changes during a hand movement task over a broad range of movement rates. We used a double readout 3D-EPI sequence at 7T to simultaneously measure VASO-CBV and BOLD responses in the sensorimotor cortex. The measured BOLD and VASO-CBV responses increased very similarly in a nonlinear fashion, plateauing for movement rates larger than 1 Hz. Our findings show a tight relationship between BOLD and VASO-CBV responses, indicating that the overall interplay of CBV and BOLD responses are similar for the assessed range of movement rates. These results suggest that the observed nonlinearity of neuronal origin is already present in VASO-CBV measurements, and consequently shows relatively unchanged BOLD responses.

Introduction
Functional Magnetic Resonance Imaging (fMRI) is the most popular means of probing neuronal activity in living humans, with blood oxygenation level-dependent (BOLD) the most common non-invasive contrast used to detect brain function. The BOLD signal relies on T2 or T2 * relaxation, which is sensitive to local concentrations of paramagnetic deoxyhemoglobin (dHb), which leads to alterations in MRI signal intensity ( Kim and Ogawa, 2012 ).
The measured BOLD signal during functional activation depends predominantly on the changes in venous blood oxygenation, which in turn depends on the induced neurovascular, hemodynamic, and metabolic changes. When interpreting functional MRI signals, we generally as- Ye ş ilyurt et al., 2008 ). In contrast, a nonlinear scaling in amplitude and duration between stimulus and the BOLD signal has been observed, for example, in visual stimuli with a short duration ( < 3-4 s). Here, BOLD response amplitude was larger than the prediction as expected from a linear system ( Birn and Bandettini, 2005 ;Liu et al., 2010 ;Logothetis et al., 2001 ;Miller et al., 2001 ;Ye ş ilyurt et al., 2008 ;Zhang et al., 2009Zhang et al., , 2008. Nonlinear behaviour was also observed for the BOLD response amplitude for a range of other fMRI experimental paradigms; word presentation rates ( Büchel et al., 1998 ;Friston et al., 1998 ;Rees et al., 1997 ), visual stimuli with different frequency contrasts ( Gaglianese et al., 2017 ). In the motor fMRI literature too, plateauing at high movement rates has been observed for different hand motor tasks ( Jäncke et al., 1998 ;Khushu et al., 2001 ;Riecker et al., 2003 ;Sadato et al., 1997 ;Siero et al., 2013 ).
Several studies combined electrocorticography (ECoG) with BOLD fMRI showing a close to linear relationship between neuronal population activity measured by electrocorticography (ECoG) and BOLD responses with respect to different stimulus rates ( Gaglianese et al., 2017 ;Siero et al., 2013 ). Specifically, in the sensorimotor cortex, Siero et al. showed that both the BOLD and neuronal ECoG responses plateau at high movement rates ( ≥ 1Hz). These findings corroborate to the suggestion that a significant part of the BOLD nonlinearity, i.e. plateauing of the response amplitude at high movement rates, has a neuronal origin ( Birn and Bandettini, 2005 ;Zhang et al., 2009 ). The remaining part is likely due to vascular nonlinear effects such as vascular refractory effects where the response to a repeated stimulus can be diminished ( Zhang et al., 2008 ). These effects are thought to occur in larger vessels that can exhibit inconsistent response delays due to pooling of upstream microvascular blood with different transit times ( Zhang et al., 2009 ). Another potential source of vascular nonlinearities is the BOLD ceiling effect, where the maximum amount of BOLD signal is reached when the blood flow is increased to the point that all the deoxyhemoglobin (dHb) in the venous vasculature is washed out ( Buxton et al., 2014 ).
Cerebral blood volume (CBV) measurements provide an alternative fMRI contrast mechanism to BOLD. The Vascular Space Occupancy (VASO) contrast is sensitive to arteriole and post-arterial CBV changes and promises higher microvascular specificity, thus better spatial localisation of the neuronal activity with reduced draining vein contamination compared to BOLD ( Hua et al., 2018 ;Huber et al., 2014 ;Jin and Kim, 2008 ;Lu et al., 2013Lu et al., , 2003. In VASO, the contrast is based on the differences between the longitudinal relaxation times (T1) of blood and tissue, and it is generated by nulling the blood signal using an inversion pulse while maintaining part of the tissue signal.
Although the linearity of the CBV response has, to date, not yet been investigated, the CBV responses are not expected to behave in the same fashion as the BOLD response to repeated stimuli and stimulus rate. Since VASO-CBV responses promise higher microvascular specificity, we might expect reduced vascular nonlinear effects as previously observed in BOLD fMRI experiments. Simultaneous assessment of CBV and BOLD responses can shed light on the extent of vascular nonlinear effects observed in the BOLD signal.
In this study, we measured simultaneously BOLD and VASO-CBV responses, using a double readout 3D-EPI VASO sequence, during the execution of hand digit movements at several movement rates. The goal was to evaluate and compare the BOLD and VASO-CBV response behaviour with respect to movement rate in sensorimotor cortex.

Subjects
Five healthy volunteers (age 25-40 years, two women and three men) without hand movement impairments participated in the study. The local ethics committee approved this study, and all volunteers provided written consent prior to participating after being informed of the experimental procedures.

Stimulus
Each subject repeated the experiment twice in two separate sessions, where the subjects performed the same motor task. The task consisted of moving the right hand periodically following a visual cue, from a rest position to a loosely clenched fist, at five different movement rates: ~0.33, 0.5, 1, 1.5, and 2 Hz. The visual cue was projected on a screen at the end of the bore of the scanner, which the subjects viewed using a mirror. The visual cue consisted of two different coloured squares (with the green square meaning move, and red indicating rest) alternating at the appropriate clenching rate and was generated in 'Matlab' (The MathWorks, Natick, United States) using the 'Psychophysics Toolbox Version 3'. Both task and movement rates have previously been shown to yield robust BOLD responses and nonlinearity with the stimulus rate ( Siero et al., 2013 ).
The task paradigm consisted of a 30 s baseline, then alternating between 12 s of movement, and 24 s of rest. In total, 12 trials were performed per movement rate in a pseudo-randomised order across subjects and sessions. All subjects were briefly trained outside the scanner and instructed to use a constant force across movement rates. The hand movements were recorded using a DataGlove 5 Ultra MRI (Fifth Dimension Technologies 5DT, www.5dt.com ) with a sampling rate of 120Hz. The right-hand dataglove is metal-free and consists of five fibre-optical sensors in total (one sensor per finger), placed on the back of each digit. These fibre-optical sensors provide an average measurement of the finger flexure for each of the five fingers by measuring the optical fibre path length. Hence, with these sensors, fist-clenching can be observed clearly, but not the movement of individual phalanges or wrist rotations. The traces from the index, middle, ring and little fingers were averaged to verify the subject's movement rate and task performance. For an example recording, see Fig. 1 . One subject's session was removed after the detection of motion-related artifacts, and the remaining 45 runs were included in the final analysis.
In order to eliminate bias in the comparison analysis between different movement rates we obtained the region of interest (ROI) from an additional stimulus run using the same paradigm described above, with 1.0 Hz and eight trials, the same processing steps were performed for all datasets including the functional localiser run.

MR sequences
Imaging was performed on a 7T MRI scanner (Philips Healthcare, Best, The Netherlands) using an 8 channel transmit coil and a 32 channel receive coil (Nova Medical Inc, Wilmington, United States) with a close to circularly polarised-mode achieved by B1-shimming of the entire brain of a group of volunteers. A Slice-Saturation Slab-Inversion (SS-SI) VASO scheme was used to simultaneously acquire CBV and BOLD weighted images using interleaved pair-wise 3D-EPI readouts ( Huber et al., 2014 ) ( Fig. 2 ). The timing parameters for the interleaved acquisition were based on previous 7T findings, taking into account gray matter (GM) and blood T 1 values, arterial arrival time in the sensorimotor cortex, and an additional margin as previously proposed ( Huber et al., 2018 ). Combining these, we used TI 1 /TI 2 /TE/TR = 1100/2600/17/3000 ms. For the magnetisation inversion, we implemented an adiabatic inversion TR-FOCI pulse, which ensures an effective inversion with reduced B1 + inhomogeneity dependency compared to a more conventional hyperbolic secant adiabatic inversion pulse ( Hurley et al., 2010 ). Data were acquired with an isotropic voxel size of 1.5 mm, FOV = 192 × 192 × 21 mm 3 , matrix size = 128 × 128, 14 slices, partial Fourier factor = 0.78 in the phase encoding direction and SENSE inplane factor = 2.5. 154 volumes were acquired per run, leading to a total scan time (per movement rate) of 7 min and 42 s.
In addition, a higher resolution T2 * -weighted scan was acquired to identify the veins in the imaging area. A 3D multishot gradient-echo EPI was acquired with the following parameters: TR/TE = 67/29 ms, flip Fig. 1. A) Single-subject dataglove traces across 12 trials shown in gray with the average response in black. The subplots illustrate the functional paradigm and task performance. From top to bottom, hand movement rates ~0.33, 0.5, 1.0, 1.5, and 2 Hz, the hand movements were performed for 12 s, followed by 24 s of rest (baseline). The task consisted of moving the right hand from a rest position (B) to a loosely clenched fist position (C). The dataglove data was used to check for task performance and to calculate an accurate movement rate for each subject. Note that, for the highest movement rates (greater than 1.0Hz), the subjects accrued timing deviations in the executed movement onset with respect to the visual cue resulting in a reduced amplitude for the average movement trace (black).

Fig. 2.
Depicted pulse sequence (sequence diagram) combined with the zmagnetisation during a 3D-EPI Slice-Selective Slab-Inversion (SS-SI) VASO acquisition. A) An adiabatic inversion (180°) pulse (TR-FOCI) followed by two acquisition modules with 3D-EPI readout with a variable flip angle ( ). Inversion times TI 1 , TI 2 , and volume repetition time TR were 1100, 2600, and 3000 ms, respectively. B) In SS-SI VASO, an increase in tissue signal is realised by manipulating the longitudinal magnetisation of the stationary gray matter (GM) tissue separately from the inflowing blood. Here, a 90°magnetisation reset pulse is applied in the imaging slab after both readout modules, resulting in an increased available tissue signal. angle = 20°, averages = 2, SENSE factor = 2, isotropic voxel size = 0.6 mm, and FOV = 215 × 160 × 75 mm 3 , matrix size = 360 × 342. 100 slices were acquired to cover the sensorimotor cortex, with a total scan time of 2 min and 35 s.

Flip angle sweep
A 3D-EPI readout not only results in higher image SNR compared a 2D-EPI readout ( Van Der Zwaag et al., 2012 ) but also leads to a reduced power deposition as typically smaller flip angles are used. In SS-SI VASO, longitudinal magnetisation recovery during the 3D-EPI readouts leads to different T1-weightings along the k z -direction, which can result in blurring across slices. To reduce this blurring, we modulated the flip angle train, such that the magnetisation in the gray matter remains approximately constant ( Gai et al., 2011 ( Fig. 2 ). The SAR values never exceeded 2 W/kg (or 20% of the local SAR limit), according to the SAR estimation of the vendor.

Data analysis
All data pre-processing was performed using the SPM12 ( Statistical Parametric Mapping) software package. We first performed a motion correction for BOLD and VASO images separately, followed by a BOLD correction used to minimise the extravascular BOLD signal contamination present in the VASO images ( Huber et al., 2014 ). No additional spatial smoothing or temporal filtering was applied to minimise loss in specificity. T2 * -weighted images were registered to the same space as the functional data, and all movement rate data sets were registered to the same common space, per subject using the functional localiser as reference.
The ROI definitions were based on VASO-CBV and BOLD activated voxels from the localiser run. After a GLM (FEAT in FSL, v.6.0) analysis, the VASO-CBV and BOLD Z-statistic activation maps were thresholded at a Z-statistic value = 2.5. To assess signals originating from draining veins, we created a "vessel ROI " based on the high-resolution T2 *weighted image, where low signal intensity spots indicate the larger draining veins, and the overlap with the BOLD activation map. To account for the extravascular BOLD signal extent of the draining vein, voxels adjacent to the large vessels were also included in the ROI ( Fig. 3 C, D in red) ( Siero et al., 2011 ). A second ROI was created based on the BOLD activation map, but excluding the overlapped voxels from the vessel ROI, resulting in a "BOLD non-Vessel ROI " (BnV ROI, Fig. 3 C, D in blue). The final ROI was defined as the common set of CBV and BOLD activated voxels, also excluding voxels from the vessel ROI, which we will dub "CBV BOLD Combined non-Vessel ROI " (CBnV ROI) ( Fig. 3 C, D in green).
To investigate the response linearity with respect to the movement rate, we first averaged the time courses of all voxels within each ROI. Next, we obtained CBV and BOLD percentage signal changes for all subjects and movement rates. To evaluate the relationship between VASO-CBV and BOLD responses in the different ROIs, we ran a linear regression on the percent signal changes per session and movement rates; five CBV-BOLD pairs of data points (per session, at five movement rates) were entered into the regression model, in total, nine slopes and intercepts were generated and further averaged. This procedure followed what was proposed in ( Makin and De Xivry, 2019 ). One-way ANOVAs with Tukey post-hoc analyses were used to assess the slope and intercept differences between ROIs. All generated plots, linear regression and ANOVA estimation, were performed using R ( R Core Team, 2020 ), where p -values < 0.05 were considered statistically significant.

Fig. 3.
Example BOLD and VASO-CBV-weighted activation maps with their respective signal amplitudes, followed by the ROI selection outline. A) shows BOLD and VASO-CBV activation maps from the functional localiser run for a representative subject, and B) signal timecourses averaged across all subjects, extracted from the ROIs used in the analysis (e.g. BnV for BOLD timecourse and CBnV for VASO timecourse). The dashed line represents the stimulation period, and the transparent shaded area represents the standard error across subjects. In C) and D), a T2 * -weighted image is shown with a schematic of the three different ROIs. In red, delineation of draining veins incorporating a BOLD extravascular signal extent (Vessel ROI). In blue, the "BOLD non-Vessel" (BnV ROI), based on BOLD activated voxels, excluding the voxels from the vessel ROI. In green, the "CBV BOLD Combined non-Vessel" (CBnV ROI), i.e. a common set of BOLD and CBV activated voxels, also excluding Vessel ROI voxels. . 1 shows the averaged dataglove response of one representative subject and illustrates the functional paradigm and associated task performance. We assessed the dataglove responses and calculated the execute movement rate for all runs. The dataglove shows that the subjects' performance was similar between the movement rates, even for higher movement rates ( > 1 Hz), where the responses start to reduce in amplitude, but the movement rate still matches with those introduced by the stimuli. Hand movements were executed within 0.337 ± 0.002 Hz, 0.505 ± 0.030 Hz, 1.002 ± 0.026 Hz, 1.523 ± 0.025 Hz, and 1.978 ± 0.114 Hz.

Fig
Robust BOLD and VASO-CBV responses in the sensorimotor cortex were detected in all participant sessions for all five different movement rates. Fig. 3 A shows exemplar BOLD and VASO-CBV activation maps for a representative subject, and Fig. 3 B the signal timecourses averaged across subjects. Fig. 4 shows the average VASO-CBV and BOLD percentage signal changes with respect to movement rate for the different ROIs. Both the BOLD and VASO-CBV response amplitude show an increase for movement rates till 1 Hz and plateaus for higher movement rates ≥ 1 Hz. As expected, the BOLD response amplitude was much higher in the vessel ROI than in the other two ROIs. VASO-CBV response amplitudes were higher in the combined ROI (CBnV ROI) than in the BnV ROI because of the way both ROIs had been defined. Interestingly, VASO-CBV responses were also significant in the vessel ROI, albeit more modest than the increase in BOLD compared to the BnV and CBnV ROIs. In addition, the VASO-CBV contrast-to-noise ratio (CNR) was approximately 56% of BOLD CNR for the CBnV ROI, calculated using the average percentage signal change across all subjects.
In Fig. 5 , a linear regression was used to assess the linearity of the relationship between VASO-CBV and BOLD responses in the three different ROIs. For all ROIs and subjects, a significant linear relationship was found between the movement rate dependent BOLD and VASO-CBV responses. The following results are the average slopes, intercepts and the Pearson's correlations across subjects and movement rates. For the vessel ROI, slope = 0.95 ± 0.10 %BOLD/%VASO-CBV, intercept = 2.20 ± 0.12 %BOLD and the R = 0.77, for the BnV ROI, slope = 1.18 ± 0.05 %BOLD/%VASO-CBV, intercept = 1.13 ± 0.07 %BOLD and the R = 0.89, and finally for the CBnV ROI, slope = 0.90 ± 0.12 %BOLD/%VASO-CBV, intercept = 1.10 ± 0.28 %BOLD and R = 0.80. The fitted linear curves are depicted in Fig. 5 for all ROIs per session; the dashed line represents the average of all nine sessions. Note that all the fitted slopes are relatively close to 1 and did not significantly differ between ROIs, as the one-way ANOVA with Tukey post-hoc analysis revealed, ROI vessel-BnV , p = 0.26, ROI vessel-CBnV , p = 0.93 and ROI BnV-CBnV , p = 0.14. Interestingly, there is a significant offset, i.e. a non-zero BOLD-intercept, at VASO CBV (%) = 0, for all ROIs. The intercept value was significantly higher for the Vessel ROI compared to both the BnV and CBnV ROIs (both p = 0.001), while the BnV and CBnV ROIs had a similar intercept value ( p = 0.9).
We also investigated the cerebrospinal fluid (CSF) contribution in the Vessel ROI. We correlated the brightness of the mean EPI from each subject, with the percentage signal change in the vessel ROI to investigate the relationship between CSF and VASO-CBV signal. The brightest voxels in the mean EPI indicates the presence of a high voxel CSF content due to the long T1 and T2 * of CSF. Lower intensity voxels generally indicate larger venous content (shorter T1 and T2 * for venous blood). A positive correlation would indicate that voxels with a high percentage of VASO-CBV signal change are related to the voxels with a large partial volume contribution of CSF. We found that the correlation for the VASO-CBV signal changes are much closer to 0 and not consistently negative or positive (i.e. non-significantly different from 0, = 0 . 44 ), suggesting that the impact of the draining veins and CSF content is on the response amplitude is small. The results for the BOLD response amplitude, however, showed a significant negative correlation ( = 7 . 71 × 10 −8 ) for all movement rates, indicating higher BOLD signal changes for voxels with higher venous blood content, as expected.

Discussion
In the present study, we evaluated the relationship between 7T BOLD and VASO-CBV fMRI responses for a simple motor task with five different hand movement rates. A double readout 3D-EPI VASO sequence enabled us to acquire BOLD and CBV responses rates in sensorimotor cortex simultaneously. As the CBV and BOLD contrast are thought to have distinct weighting for (micro)vascular compartments, we anticipated distinct response behaviour with respect to movement rate. Interestingly, we found a high degree of similarity between BOLD and CBV response behaviour regarding movement rate.
Observed VASO-CBV activation patterns ( Fig. 3 A) are similar in magnitude and spatial extent to previous studies using the same spatial resolution in sensorimotor cortex, showing reduced large draining vessel activations compared to the BOLD activation patterns ( Huber et al., 2018( Huber et al., , 2014. The relative sensitivity differences between VASO-CBV and BOLD was also in agreement with previous work ( Huber et al., 2014 ). For the CBV BOLD Combined non-Vessel ROI (CBnV ROI), the VASO CBV contrast-to-noise ratio (CNR) was approximately 56% of BOLD CNR.
We found that the BOLD response first increases and then plateaus for movement rates ≥ 1 Hz. Interestingly, the CBV response showed a markedly similar response behaviour as the BOLD response for all ROIs ( Fig. 4 ). Several fMRI studies have described movement frequency dependency in sensorimotor areas and other areas such as basal ganglia, cerebellum and spinal cord ( Agnew et al., 2004 ;Jäncke et al., 1998 ;Fig. 4. The BOLD and VASO-CBV percent signal change with respect to movement rate for the different ROIs, averaged across subjects. A.) Vessel ROI, B.) The "BOLD non-Vessel " (BnV ROI), and C.) The "CBV BOLD Combined non-Vessel " (CBnV ROI, i.e. common CBV and BOLD activated voxels). The higher sensitivity of BOLD compared with VASO-CBV is reflected in the higher percentage signal change. The error bars are the standard error. Regardless of the ROI, both BOLD and VASO-CBV percentage signal change show similar behaviour, with both increasing in a linear fashion for movement rates < 1 Hz and saturating for movement rates ≥ 1 Hz.

Fig. 5. %BOLD versus %VASO-CBV for different ROIs, showing the linear trend in the relationship between both measured responses for all movement rates.
Each line represents the slope of %BOLD versus %VASO-CBV for all movement rates (open circles) of a single session, and the dashed lines represent the average slope across all sessions. In red, the draining vein or "vessel" ROI. In blue, the "BOLD non-Vessel" (BnV) ROI; based on BOLD activated voxels, excluding the voxels contained in the vessel ROI. In green, the "CBV BOLD Combined non-Vessel" ROI (CBnV ROI), i.e. the common set of BOLD and VASO-CBV activated voxels excluding vessels. Maieron et al., 2007 ;Rao et al., 1996 ). In addition, electrophysiology studies have shown evidence for a negative relationship between motor cortex neuronal firing rate and movement frequency: greater speed of movement was associated with a reduced number of spikes/second in primates ( Aflalo and Graziano, 2007 ). This notion is in agreement with human electrophysiology (ECoG) studies which observed that the plateauing could be associated with strongly reduced neuronal population activity (high-frequency gamma-band power) for fast movement rates after the very first movement. Moreover, the beta-band power re-mained suppressed, indicating a release of motor cortex inhibition, in a sense facilitating movement initiation for these fast movement rates ( Hermes et al., 2012 ;Siero et al., 2013 ).
Plateauing of the BOLD response has been observed in sensorimotor cortex using a range of task designs such as finger tapping, buttonpressing or hand movements ( Jäncke et al., 1998 ;Khushu et al., 2001 ;Riecker et al., 2003 ;Sadato et al., 1997 ;Siero et al., 2013 ). In all cases, response plateauing was observed at high movement rates with a variation in the exact range of frequencies, differences that are likely related to the different experimental designs.
Interestingly, the observed VASO-CBV signal changes in the vessel ROI were comparable to the VASO-CBV signal changes in the CBnV ROI and higher than the BnV ROI. A possible explanation could be a partial volume effect from cerebrospinal fluid (CSF) ( Donahue et al., 2006 ;Huber et al., 2015 ). The CSF signal is much higher than the gray matter signal for the used SS-SI-VASO scan parameters. Thus, for a similar displacement of the voxel's content by vascular dilation and/or adjacent tissue swelling, the signal intensity is reduced much more for high CSF content voxel than for low CSF content voxel. We investigated the CSF contamination by performing a voxelwise correlation between the mean EPI signal intensity and %VASO-CBV and %BOLD signal changes. While we found a negative correlation in the BOLD data confirming the importance of veins in this ROI, no such correlation was found for the VASO-CBV contrast, indicating minimal CSF contamination in the vessel ROI. The result of the correlation analysis is shown in the supplementary material (Suppl. Fig. 1). However, other potential confounders could also explain the unexpected vessel ROI VASO CBV signal change. Previous studies reported that VASO fMRI is expected to capture CBV changes in small pial arteries ( Donahue et al., 2006 ;Huber et al., 2015 ). Small pial arteries are most likely present in the vessel ROI but exhibit no image contrast, as the pial veins do, with respect to the CSF and gray matter signal. We also checked the VASO-CBV percentage signal change in the vessel ROI and the CBnV ROI when omitting the BOLD contamination correction. The results are shown in the supplementary material (Suppl. Fig. 2). As expected, the VASO-CBV responses are reduced when omitting the BOLD contamination correction. The amplitude for the non-BOLD corrected vessel ROI responses are comparable to the non-BOLD corrected CBnV ROI responses for all frequencies besides 1.5Hz, similar to the corrected data. Hence, any BOLD contamination is not likely to be the cause of the comparable response amplitude for the vessel ROI with respect to the CBnV ROI.
The observed tight relationship between BOLD and VASO-CBV responses for all ROIs indicates that blood volume changes are strongly associated to changes in deoxyhemoglobin, at least over the range of examined movement rates. This can also be seen in the scatter plot and fitted linear equations depicted in Fig. 5 , showing a high correlation between %BOLD and %VASO-CBV changes and an almost one-to-one relationship as seen from the fitted slopes for all ROIs. Interestingly, we observed a non-zero intercept for the fitted linear equations describing the %BOLD versus %VASO-CBV all ROIs. At first instance, a zerointercept would be expected, meaning that in the absence of a notable CBV change, a zero % BOLD signal change would be observed. The observed non-zero BOLD intercept, though, could be remaining extravascular BOLD signal in the VASO-CBV data, i.e. residual BOLD contamination. Both BOLD and VASO-CBV readouts have similar T2 * weighting, in this case, the positive BOLD signal changes cancel some of the negative VASO-CBV signal changes during activation. Most of this T2 * or BOLD contamination is expected to be removed from the generated VASO-CBV images using a BOLD correction scheme ( Huber et al., 2014 ); however, remaining contamination cannot be fully ruled out. It is worth pointing out that the BOLD correction scheme has been validated by at least two different designs, in both cases, correcting major part of the BOLD contamination with minimal overestimation of VASO-CBV signal changes ( Huber, 2015 ). Another possibility is that the BOLD intercept reflects an (intracortical) draining vein effect, which is presumably not present in the VASO-CBV data but likely is present in the BOLD data, even after masking out the largest vessels visible on the high-resolution T2 * anatomy image. Note that the intercept was the highest for the Vessel ROI.
The nonlinear behaviour observed in both the BOLD and VASO-CBV responses regarding movement rate is most likely, and for the most part, explained by the neuronal source ( Hermes et al., 2012 ;Siero et al., 2013 ). The remaining part is probably of vascular origin and caused by vascular nonlinearity effects such as vascular refractory effects and the ceiling phenomenon. The ceiling effect could potentially be the cause of the BOLD nonlinearity, as suggested previously ( Birn and Bandettini, 2005 ;Bruhn et al., 1994 ;Buxton et al., 2004 ). However, other fMRI modalities such as CBF mapping using Arterial Spin Labelling (ASL) have found a comparable nonlinear response with respect to different stimulus duration, and a linear relationship between CBF and BOLD ( Miller et al., 2001 ), which are also similar to previous PET studies ( Blinkenberg et al., 1996 ;Dettmers et al., 1996 ;Sadato et al., 1997 ). Our CBV findings complement these previous findings, and together they suggest that there is a minimal contribution of the ceiling effect in the BOLD nonlinearity since BOLD response is driven by CBF and CBV changes ( Buxton et al., 2014 ). A plausible explanation for the VASO-CBV nonlinearities is that the CBV measurements reflect the nonlinear transformation from the stimulus to the neural response since CBV are closely related to the underlying neuronal activity. Consequently, the BOLD nonlinearities could be an extension of that, as the CBV nonlinearities are also part of the BOLD responses.
We observed a lower sensitivity for VASO-CBV responses, VASO-CBV CNR was about 56% compared to the BOLD CNR. The higher sensitivity is expected for GE-BOLD due to the nature of the contrast mechanism. In addition, residual signals from macrovessels or signals from non-specific large draining veins, can increase the observed BOLD CNR ( Huber et al., 2017 ;Uluda ğ and Blinder, 2016 ). With respect to the sequence used here, the SS-SI 3D-EPI sequence provides simultaneous acquisition of VASO-CBV and BOLD fMRI measurements, however, care should be taken regarding the comparison and interpretation. An optimal VASO sequence requires careful adjustment of timing parameters; an appropriate inversion time that is not too long to avoid inflow effects, and a long TR to satisfy the VASO blood refilling conditions; usually TR = 3 s or longer to allow blood to refill the vasculature in the imaging volume but also to allow leaving the imaging volume before the next inversion pulse ( Huber et al., 2014 ;Jin and Kim, 2008 ). In addition, a short TE is required to minimise BOLD contamination. In contrast, an optimal GE-BOLD sequence requires a relatively longer TE, and the TR is usually limited by the coverage of the field of view. Thus, the current TR restriction of the optimal SS-SI 3D-EPI VASO limits a more thorough investigation of the temporal features of the hemodynamic response and the neurovascular coupling. Methods like Multiple Acquisitions with Global Inversion Cycling (MAGIC) ( Lu et al., 2004 ;Scouten and Constable, 2007 ) and the more recently Multiple Acquisitions with Global Excitation Cycling (MAGEC) ( Huber et al., 2020 ) can achieve whole-brain coverage and, therefore, be used to increase the temporal resolution in future VASO studies.

Conclusion
We observed a strong linear relationship between VASO-CBV and BOLD responses, both similarly increasing with increasing hand movement rate and plateauing at high movement rates ≥ 1 Hz. The presumed higher microvascular specificity of VASO-CBV compared to BOLD does not directly result in a more linear response behaviour at high hand movement rates. The observed plateau behaviour, i.e. nonlinear effects, regarding movement rate, are likely predominantly of neuronal origin and therefore present in both the VASO-CBV and BOLD response.

Data availability
The data that support the findings of this study are available from the corresponding author, upon reasonable request.