Breath-hold BOLD fMRI without CO 2 sampling enables estimation of venous cerebral blood volume: potential use in normalization of stimulus-evoked BOLD fMRI data

BOLD fMRI signal has been used in conjunction with vasodilatory stimulation as a marker of cerebrovascular reactivity (CVR): the relative change in cerebral blood flow (CBF) arising from a unit change in the vasodilatory stimulus. Using numerical simulations, we demonstrate that the variability in the relative BOLD signal change induced

end-tidal carbon dioxide, unlike traditional methods to estimate BOLD-CVR.It also more closely represents a specific physiological characteristic of brain vasculature than BOLD-CVR, namely blood volume.

Introduction
The signal measured in a blood oxygenation level dependent (BOLD) contrast functional magnetic resonance imaging (fMRI) experiment largely reflects the change in the paramagnetic deoxygenated hemoglobin content of the image voxels.It has been shown that two modulators of the BOLD signal explain much of its variance between subjects: the cerebrovascular reactivity (CVR) and baseline deoxyhemoglobin content (Liu et al., 2013).CVR is normally considered as the relative increase in the cerebral blood flow (CBF) in response to a vasodilatory stimulus, where higher CVR leads to increased BOLD signal changes throughout the brain of healthy subjects (Handwerker et al., 2007;Liau and Liu, 2009;Liu et al., 2013;Thomason et al., 2007).Baseline deoxyhemoglobin content, instead, depends on the oxygen extraction fraction and the blood volume responsible for generating the BOLD signal (CBV dHb ) and is directly proportional to the BOLD signal change (Buxton et al., 2004;Davis et al., 1998;Hoge et al., 1999;Liu et al., 2013;Lu et al., 2010Lu et al., , 2008)).The deoxygenated blood is found principally in the capillaries and veins, the volume of which, in a given image voxel, depends not only on the vascularity of the local tissue but also on the fraction of the voxel volume occupied by brain tissue.This so-called partial volume varies considerably when voxel volumes are large, as they typically are for fMRI studies.The dependence of the BOLD signal response to these underlying physiological factors may be exploited to characterize brain tissue physiology.
CVR can be measured directly by applying a vasodilatory stimulus that increases CBF (Hall et al., 2014;Wei et al., 1980).Common vasodilatory stimuli include acetazolamide (Vagal et al., 2009) or hypercapnia (increase in arterial carbon dioxide (CaCO 2 )) induced by breathing exogenous carbon dioxide (CO 2 ) or induced endogenously through breath-holding (Bright and Murphy, 2013;Kastrup et al., 2001;Moreton et al., 2016).Parametric maps of CBF, and thus CVR, can then be calculated based on dynamic susceptibility contrast or arterial spin labeling MRI (Nöth et al., 2006;Tancredi and Hoge, 2013).CVR has also been investigated by exploiting the dependence of the BOLD signal on relative changes in CBF in response to, what is assumed to be, isometabolic hypercapnia (Bulte et al., 2012;Davis et al., 1998;Hoge et al., 1999;Liu et al., 2019), here termed BOLD-CVR: Here, TE is the echo time of the imaging sequence, A and β are constants that depend on the vessel size and geometry and the strength of the magnetic field, α is the Grubb exponent (Grubb et al., 1974), CBV dHb is the non-arterial cerebral blood volume (approximately 70% of the total CBV) responsible for generating the BOLD signal, SvO 2 is the venous oxygen saturation, and [Hb] is the hemoglobin concentration in blood.Zero subscripts denote baseline values.Eq. ( 1) can be linearized with respect to relative changes in CBF and reformulated as a function of CVR as (see the Appendix): where ΔCaCO 2 denotes a change in arterial carbon dioxide.Eq. ( 2) highlights that changes in the BOLD signal depend on the CVR and the voxel baseline deoxyhemoglobin content, itself determined by the CBV, oxygen extraction fraction (OEF), and blood hemoglobin content.Parametric maps of BOLD-CVR are calculated by normalizing the BOLD signal modulation in Eq. ( 2) by the change in the end-tidal CO 2 (PetCO 2 ).Breath-holding offers a convenient alternative to administering CO 2 during MRI (Bright and Murphy, 2013;Kastrup et al., 2001;Moreton et al., 2016).Consequently, breath-hold-induced variations in the BOLD fMRI signal have been used as a proxy of CVR.During breath-holding, regional fractional changes in CBF are relatively uniform across the brain and they are reproducible across subjects (Taneja et al., 2020).It follows, therefore, from Eq. ( 2) that a breath-hold fMRI experiment may enable us to estimate a marker of CBV dHb (the non-arterial blood volume), which we term BOLD-CBV to highlight its derivation from BOLD measurements.Such maps could be obtained by normalizing the breath-hold BOLD signal modulation by the purely intravenous BOLD signal, for example, from a large vein where CBV dHb approaches unity.This approach would not require measurements of PetCO 2 , thus eliminating a potential source of inaccuracy.Notably, the normalization by a venous signal, differently from the normalization by the change in PetCO 2 , accounts for all other physiological parameters except that of interest, assuming that these parameters are reasonably uniform across the brain.In particular, the observation of a uniform CVR across the brain is supported by arterial spin labeling-based measurements (Taneja et al., 2020).Finally, the spatial variability of CBV dHb across the brain is expected to be higher than that of CVR (Zhao et al., 2021).
The physiological factors affecting the BOLD signal response vary between individuals and across the brain.As such, they may be considered a nuisance because of the variability they add to BOLD response data in neuroscientific fMRI studies, in which it is desirable to maximize the sensitivity to detecting task-related changes in brain activity in a cohort of participants.Therefore, efforts have been made to minimize the effects of these unwanted sources of inter-subject variability by measuring and accounting for them, thus increasing the sensitivity and specificity for detecting neural activity-related BOLD signal changes.Previous inter-subject normalization approaches have employed BOLD-CVR or CBF measurements, hypercapnic challenges, resting state fMRI experiments, or the residuals of task-evoked fMRI (Bandettini and Wong, 1997;Cohen et al., 2004;Kalcher et al., 2013;Kannurpatti and Biswal, 2008;Kazan et al., 2016;Murphy et al., 2011).Notably, BOLD-CVR is a practical approach for correcting vascular sources of between-subject differences in the stimulus-evoked BOLD response, leading to more consistent activation patterns across participants in task-relevant regions (Murphy et al., 2011).However, estimating CVR through the BOLD signal (BOLD-CVR) suffers from two problems: (1) it is confounded by the BOLD signal dependence on the underlying physiology, in particular, CBV dHb , but also OEF and [Hb], and (2) estimating CVR requires reliably measuring the extent of the vasodilatory stimulus which, for hypercapnia, necessitates accurately measuring PetCO 2 .
Firstly, the present study aimed to develop and test an approach for calculating BOLD-CBV, a marker of non-arterial blood volume (denoted in Eqs. ( 1) and (2) as CBV dHb ).We identified the dependency of the breath-hold-induced BOLD signal change on CBV dHb allowing the calculation of a physiological parameter that did not require measurement of PetCO 2 .Using numerical simulations, the theoretical dependency of BOLD-CBV on other vascular parameters influencing the BOLD signal change was investigated.
As a potential application of the BOLD-CBV marker, the present study secondly compared the extent to which voxel-wise measures of vascular characteristics, BOLD-CVR and BOLD-CBV, accounted for between-subject variability in the BOLD response to a visuomotor task in a higher-level fMRI analysis.

Simulations
Simulations aimed to investigate the relationship between the BOLD fMRI signal during a breath-hold and each of the following vascular parameters: the blood volume responsible for generating the BOLD signal at baseline (CBV dHb,0 ) (Chiarelli et al., 2022a), the relative change in CBF during breath-holding (that is, the CVR), OEF 0 , and [Hb].In addition, we aimed to establish a practical marker of deoxygenated blood volume based on normalizing the breath-hold induced tissue BOLD fMRI signal by the signal arising from the intravascular venous compartment such as the superior sagittal sinus (SSS).Such a normalization approach was similar to one previously proposed to calculate CBV based on hyperoxic BOLD fMRI data (Bulte et al., 2007).We term this blood volume marker BOLD-CBV.Our simulations also compared BOLD-CBV with CBV dHb to assess if BOLD-CBV was a good marker of CBV dHb .This local marker of blood volume, BOLD-CBV, assumes a uniform fractional increase in CBF during breath-hold, namely that the change in OEF in the SSS closely represents that across all brain tissue (Zhao et al., 2021).
The BOLD signal change resulting from hypercapnic vasodilation was simulated using a multi-compartment vascular model in which Monte Carlo simulations of the apparent transverse relaxation rate (R 2 *) effects were combined with empirical observations of intravascular R 2 *, both due to the presence of deoxyhemoglobin in blood, and estimated at 3 T.One extravascular and three intravascular compartments were considered.Vascular contributions comprised the arteriole, capillary, and venule compartments, all occupying a given percentage of the imaging volume (CBV a , CBV c , and CBV v , respectively) (Uludag et al., 2009).All blood volumes were expressed as a function of CBV dHb (Eq.(3)), defined as the sum of the capillary and venous blood volumes, CBV dHb = CBV c + CBV v .In the model, a random variable ρ was defined, representing the ratio of CBV dHb to CBV c .CBV a was empirically set equal to 50% of CBV c , drawing from previous work and noting that the arteriolar contribution to the BOLD signal is small (Uludag et al., 2009).Hence, the total blood volume CBV tot was equal to: The Object-oriented Development Interface for Nuclear Magnetic Resonance environment (Jochimsen and von Mengershausen, 2004) was used to simulate extravascular R 2 * effects due to deoxyhemoglobin.In detail, the simulations added the effects generated from vessels modeled as cylinders with two radii (4 and 20 μm), and the susceptibility frequency shifts at the surface induced by the variable amount of deoxyhemoglobin in the vessels.
Random vessel orientations were accounted for by varying the angle between the main field and the cylinder between 0 and π/2 rad.The small and large vessel radii were respectively chosen to reflect in vivo reports of capillary radius (4 μm) and the mean of the arteriolar (17.8 μm) and venular radii (23.3 μm), which was equal to 20 μm (Adams et al., 2015;Saiga et al., 2021).
The frequency shift at the vessel surface was calculated as Uludag et al. (2009): where Hct [%] = 3*Hb [g/dl] is the hematocrit, SO 2,off is the hemoglobin oxygen saturation that produces no shift between intravascular and extravascular fluids, SaO 2 , SvO 2 or ScO 2 are the hemoglobin oxygen saturations of the arterial, venous or capillary vessels, B 0 is the static magnetic field, Δχ do is the magnetic susceptibility shift between fully oxygenated and fully deoxygenated red blood cells, and γ is the proton gyromagnetic ratio.The extravascular compartment also considered water diffusion with a diffusion coefficient of D = 1 μm 2 /ms.
The deoxyhemoglobin-induced extravascular effects on R 2 * were added onto the baseline tissue apparent R 2 * (in Hz), which was assumed to be Uludag et al. (2009): hence, the total extravascular signal was estimated as: In each vascular compartment (a = arterioles, c = capillaries, v = venules), intravascular relaxation rates were calculated as Zhao et al. (2007): with the total intravascular signal being: The total BOLD signal was assumed to be: The arterial oxygen saturation SaO 2 was calculated as a function of the oxygen partial pressure (PaO 2 ) according to the Severinghaus equation: Based on the definition of the OEF (Appendix, Eq. ( A1)) and based on the equation for oxygen concentration in blood (Appendix, Eq. (A2)) the following expression for venous saturation is obtained: The capillary oxygen saturation (ScO 2 ) was instead estimated as: Assuming an isometabolic vasodilation during breath-hold and based on Fick's principle, the OEF during breath-hold was inferred as: The increase in CBV was estimated using the Grubb relation linking fractional changes in CBV with fractional changes in CBF: with α being the Grubb exponent.
Finally, the relative BOLD signal change in response to breath-hold was set equal to: The physical and physiological parameters of the simulation model are shown in Table 1 and were either fixed or considered as random variables with their statistical distributions chosen to simulate the variability found in vivo (Chiarelli et al., 2022a).In particular, the range of CBV dHb values was set in line with Eq. (3) (the relationship between CBV dHb and the total CBV) based on previously reported ranges of CBV in healthy subjects (Leenders et al., 1990).The range of CBF/CBF 0 values was chosen in accordance with previously reported ranges of CVR in healthy subjects (Zhao et al., 2021).The difference between the arterial oxygen tension at baseline (PaO 2,0 ) and breath-hold (PaO 2 ) was chosen to reflect the average trace of end-tidal O 2 measured across subjects in vivo, which varied between 108 and 123 mmHg over time.Modeling the difference between PaO 2,0 and PaO 2 enabled us to account for the mild hypoxemia induced by a breath-holding task.Of note, different from Eq. (2) in the Introduction, which was a simplified representation of the dependency of the BOLD signal on various physiological parameters, complete arterial oxygen saturation was not assumed in the numerical simulations, as PaO 2 in Eq. ( 10) was modeled as a random variable.
Correlations and least-squares fits were then calculated between each vascular parameter, and both the BOLD signal and the BOLD signal normalized to its intravascular component.

Participants and MRI data acquisition
Twenty-seven right-handed healthy subjects (38.1 ± 11.0 years old, 15/12 females/males) underwent MRI on a GE Signa HDx 3 T system (GE Medical Systems, WI, US) using an 8-channel receive-only head coil.All subjects gave written informed consent.The imaging protocol included structural T 1 -weighted MRI, stimulus-evoked BOLD fMRI, breath-hold BOLD fMRI, and B 0 inhomogeneity field map acquisitions with the parameters shown in Table 2.
During stimulus-evoked BOLD fMRI, subjects performed an upper limb visuomotor serial reaction time (SRT) task, which requires the association and integration of several domains underlying motor sequence learning (Nissen and Bullemer, 1987).Here, the position of a visual stimulus was reported by the participant by pressing the corresponding key on a keypad.This same protocol was previously applied in a cohort of patients (Lipp et al., 2021).Briefly, the task comprised a series of Sequence, Random, and Rest blocks.Sequence blocks included four repeats of eight stimuli and lasted 16 s.Random blocks presented stimuli in a random order, with no stimulus appearing twice consecutively, and lasted 16 s.One Random block was presented every three Sequence blocks.Rest blocks presented a static image matching the task and lasted on average 15 s.A visual representation of the SRT task is shown in Supplementary Fig. 1.Overall, the experiment took 7 min and 6 s to complete.
For breath-hold BOLD fMRI, the breath-hold task was adapted from Murphy et al. (2011).Continuous CO 2 partial pressure signals in the expired air were measured throughout the experiment via a nasal cannula (AEI Technologies, PA, US).Following an initial baseline of approximately 30 s during which no task was performed, subjects were given instructions on a screen and guided through four cycles of breath-holding and recovery.Each cycle comprised the following phases: 18 s of paced breathing (alternating breathing in and breathing out for 3 s each), 15 s of end-expiration breath-holding, exhalation, and 15 s of recovery (spontaneous breathing with no specific instructions).Overall, the experiment took 3 min and 45 s to complete.

PetCO 2 calculation
For each subject, a trace of PetCO 2 was calculated and considered a surrogate measure of arterial CO 2 changes (McSwain et al., 2010;

Table 1
Variables used in the simulations are reported in alphabetical order.N and Γ respectively denote the normal and gamma distributions for random variables.Takano et al., 2003).In particular, the PetCO 2 trace was calculated by isolating and linearly interpolating the peaks of the raw CO 2 trace (Fig. 1) (Chiarelli et al., 2022b).

Image analysis 2.4.1. Breath-hold fMRI: BOLD-CVR and BOLD-CBV
The processing steps illustrated in this section are summarized in Fig. 1.
For each subject, breath-hold BOLD fMRI images were processed using FMRIB's linear image registration tool for motion correction (MCFLIRT) (Jenkinson et al., 2002), expressed voxel-wise as fractional BOLD changes relative to their temporal mean and high-pass Gaussian filtered (full-width half maximum = 60 s).
To enable voxel-wise estimation of BOLD-CVR, the PetCO 2 trace was resampled to match the number of time points in the breath-hold fMRI acquisition and detrended via quadratic polynomial fitting (Chiarelli et al., 2022b).The resampled and detrended PetCO 2 trace was shifted in time via convolution with a canonical hemodynamic response function, which was generated using spm_hrf from the Statistical Parametric Mapping toolbox (v12) (Friston et al., 2007).A map of BOLD-CVR expressed as the BOLD signal percent change per unit of PetCO 2 was then calculated using the General Linear Model framework (FSL film_gls) and the subject-wise PetCO 2 trace as a regressor (Chiarelli et al., 2022b).
To enable voxel-wise estimation of BOLD-CBV, a region of interest (ROI) was manually delineated using ITK-SNAP (Yushkevich et al., 2006) in the posterior SSS based on four contiguous sagittal slices of the T 1 -weighted structural image.A map of gray matter (GM) partial volume was obtained from the T 1 -weighted structural image using FMRIB's Automated Segmentation Tool (FAST) (Zhang et al., 2001).The breath-hold BOLD fMRI image was aligned with the T 1 -weighted image using FMRIB's linear image registration tool (FLIRT) (Jenkinson et al., 2002;Jenkinson and Smith, 2001).The resulting transformation matrix was inverted to bring the T 1 -weighted image, the map of GM partial volume (both using trilinear interpolation), and the SSS ROI (using nearest-neighbor interpolation) into breath-hold BOLD fMRI space.The GM partial volume map in breath-hold BOLD fMRI space was thresholded at 0.5 and binarized to obtain a GM segmentation.
In breath-hold BOLD fMRI space, an estimate of the breath-hold BOLD fMRI signal time-course in the SSS (SSS BH-BOLD ) was obtained as follows, aiming to establish the intravascular BOLD signal variation while limiting the influence of partial volume effects on this estimate.Within the co-registered SSS mask, the signals having a high covariance (≥ 90th percentile) with the mean GM breath-hold BOLD fMRI signal and having a baseline BOLD signal value between 20% and 100% of the average GM baseline BOLD signal were averaged.The voxels that covaried strongly with the global GM signal were selected to eliminate voxels with limited relevant BOLD variability, probably caused by partial volume effects within the large voxels of the functional images.
Including voxels with a BOLD signal baseline intensity within a plausible range for vascular tissue eliminated large fractional oscillations in the signal, which were probably induced by the effects of motion at the border of the SSS.BOLD-CBV (fractional BOLD-generating blood volume) was calculated voxelwise by regressing the breath-hold BOLD fMRI signal against SSS BH-BOLD .Hence, the resulting BOLD-CBV estimate was expressed in quantitative units and could be compared across subjects.The amplitude of BOLD signal oscillations synchronous to the breath-hold task was also evaluated by multiplying the BOLD-CBV by the standard deviation of the breath-hold signal in the SSS.
BOLD-CBV and BOLD-CVR maps calculated in one representative subject are shown in Fig. 2A-B.

First-level analysis of stimulus-evoked fMRI
For each subject, brain extraction was performed on the T 1 -weighted structural image using the FSL brain extraction tool (BET) (Smith, 2002) with bias field correction and neck clean-up and a fractional intensity threshold of 0.2.These parameters were chosen to suit this cohort of subjects.
A B 0 field map was calculated as the complex division between the gradient-recalled echo (GRE) echo planar imaging (EPI) data acquired at two TEs.The first-echo magnitude image was calculated as the absolute value of the complex GRE image at TE = 7 ms.Brain extraction on the first-echo magnitude image was performed using BET with robust brain center estimation and a fractional intensity threshold of 0.5.
Brain activation in response to the SRT task was evaluated by processing the stimulus-evoked BOLD fMRI images using the FSL fMRI expert analysis tool (FEAT) and a first-level analysis (Woolrich et al., Fig. 1.Processing pipelines for BOLD-CVR and BOLD-CBV.This scheme illustrates the processing steps employed to calculate: (1) a map of BOLD-CVR based on the CO 2 signal and breath-hold BOLD data and (2) a map of BOLD-CBV based on the breath-hold BOLD and T1w data.Abbreviations.3D/4D: three/four-dimensional; BH: breath-hold; BOLD: blood oxygenation level-dependent; CBV: cerebral blood volume; CO 2 : carbon dioxide; CVR: cerebrovascular reactivity; GM: gray matter; HRF: hemodynamic response function; PetCO 2 : partial pressure of end-tidal carbon dioxide; SSS: superior sagittal sinus; T1w: T 1 -weighted; ROI: region of interest.
E. Biondetti et al. 2004Biondetti et al. , 2001)).B 0 unwarping was performed based on the B 0 field map (Jenkinson, 2004(Jenkinson, , 2003)), followed by registration to the T 1 -weighted structural image using FLIRT and registration from the T 1 -weighted structural image to the 2-mm isotropic Montréal Neurological Institute (MNI) standard space using FLIRT and FMRIB's nonlinear image registration tool (FNIRT) (Andersson et al., 2007a(Andersson et al., , 2007b)).MCFLIRT was applied to functional data, followed by non-brain tissue removal using BET, spatial smoothing using a Gaussian kernel (full width at half maximum = 5 mm), grand-mean intensity normalization by a single multiplicative factor, and high-pass temporal filtering using a Gaussian-weighted least-squares line (sigma = 50 s).Time-series statistical analysis was performed using FMRIB's improved linear model (FILM) with local autocorrelation correction.The stimulus-evoked BOLD fMRI signal was modeled with a block task regressor plus its temporal derivative (Friston et al., 1998;Woolrich et al., 2001).

Higher-level analysis of stimulus-evoked fMRI
To ensure alignment between the breath-hold and stimulus-evoked BOLD fMRI data in MNI space, the BOLD fMRI data were distortioncorrected based on the B 0 field map.The same correction was applied to the maps of BOLD-CVR and BOLD-CBV.These maps were then transformed to the 2-mm isotropic MNI standard space by applying a transformation obtained by concatenating a FLIRT transformation from the distortion-corrected breath-hold BOLD fMRI to T 1 -weighted space and a FLIRT plus FNIRT transformation from T 1 -weighted to MNI space.BOLD-CBV and BOLD-CVR maps averaged across subjects in MNI space are shown in Fig. 2C-D.
To empirically assess which vascular covariate, whether BOLD-CVR or BOLD-CBV, best represented the main source of inter-subject variability in BOLD fMRI, higher-level FEAT analyses were performed using three models: the first modeled the average task response only (SRT); the second and third added a vascular explanatory variable modeling the voxelwise variation in BOLD-CVR (SRT+BOLD-CVR) or BOLD-CBV across subjects (SRT+BOLD-CBV), respectively.
To enable a regional analysis of task-associated brain activation in each higher-level model, four ROIs were selected, encompassing regions involved in SRT-related motor learning and processing of visual stimuli (Hardwick et al., 2013).The ROIs included: the precentral gyrus, which contains Brodmann's area 4 and corresponds to the primary motor cortex; the juxtapositional lobule (or supplementary motor area), which contains Brodmann's area 6 and is involved in motor planning; the inferior division of the lateral occipital cortex, which contains Brodmann's areas 18 and 19 and corresponds to the secondary (V2) and associative visual cortices (V3-V5); and the superior parietal lobule, containing Brodmann areas 5 and 7, which are respectively responsible for visuomotor coordination and somatosensory association.ROI delineation (Fig. 3) was based on the 2 mm-isotropic Harvard-Oxford cortical atlas thresholded at a probability equal to 0.25, which was observed to produce no overlap between contiguous regions.
To evaluate the effect of each higher-level model on the extent of task-associated brain activation in each ROI, a histogram analysis was performed.The following were assessed: (1) the Z score distribution corresponding to the task response for each of the three models; (2) the Z score distribution corresponding to the association between task response and vascular covariate in the SRT+BOLD-CVR and SRT+BOLD-CBV models.MATLAB's histfit.m function was used, setting the number of bins equal to the square root of the number of elements in the data set and fitting a nonparametric kernel-smoothing distribution to the histogram.For each model and each ROI, this analysis yielded: (1) the strength of activation in response to the task, and (2) the number of voxels where an association between the task response and one of the two vascular covariates was found.Both ( 1) and ( 2) were evaluated based on the number of voxels in the ROI taking on a range of significant Z score values (i.e., Z > 3.1).For both ( 1) and ( 2), the frequency distribution of the number of voxels in each ROI was displayed using histogram plots.
To assess the generalizability of the higher-level analysis, bootstrapping was performed.In brief, bootstrapping enables inference about a population based on sample data by randomly resampling the sample data with repetition multiple times, each time deriving the same metrics from the resampled population.A distribution of values can then be derived for each metric, based on which inference can be performed.By bootstrapping with 1000 iterations, for each model (SRT, SRT+BOLD-CVR, SRT+BOLD-CBV) an estimate of the distribution was calculated for (1) the average number and Student's T statistic of significantly activated voxels in response to the SRT task, (2) the average number and Student's T statistic of voxels in which the association between the response to the task and a vascular covariate was significant.Differences between the three models were tested by calculating pairwise differences between the values in (1) and (2).

Data and code availability statement
In line with local ethics guidelines and subject privacy policies, the data and code supporting the findings in this study are available via request to the authors.Data are subject to a formal data-sharing agreement.

Variables affecting the breath-hold BOLD signal
Fig. 4 shows the Monte Carlo simulation results.The BOLD signal and the BOLD signal normalized to its intravascular component (i.e., BOLD-CBV) were strongly and positively correlated with CBV dHb (Fig. 4A, E).The BOLD signal was also positively, although less strongly, correlated with OEF 0 (Fig. 4B), the hemoglobin level (Fig. 4C), and the relative change in CBF (Fig. 4D), representing CVR.In the case of hemoglobin, the correlation was weak (Fig. 4C).The BOLD-CBV signal was negatively correlated with OEF 0 (Fig. 4F) and positively, although weakly, correlated with both the hemoglobin concentration (Fig. 4G) and the relative change in CBF (Fig. 4H).For data acquired in vivo, in the GM the average BOLD-CBV was positively and significantly correlated with the average BOLD signal changes induced by the breath-hold task (r = 0.42, p = 0.03).This result was in line with the simulations shown in Fig. 4A.

Inclusion of vascular covariates in group-level fMRI
Fig. 6 shows the group response to the SRT task estimated by the three models (A-L) and the association between the task response and each vascular covariate (M-T).Fig. 7 reports the group-level statistical distribution of Z scores corresponding to the association between the task response and each vascular covariate in each ROI.The statistical distribution of Z scores corresponding to brain activation in response to the task in each ROI is shown in Supplementary Fig. S2.
Both group-level corrections for vascular covariates SRT+BOLD-CVR and SRT+BOLD-CBV moderately increased the significance of activation to the task, especially in the motor cortex at superior levels (Fig. 6B-D, 6F-H, 6J-L, and 7A-B).The largest increase was observed in SRT+BOLD-CBV (Figs. 6 and 2).Moreover, the BOLD-CBV covariate showed a spatially larger and stronger correlation with the task response, especially in the precentral gyrus and superior parietal lobule (Figs.6F-H,  7A-B, and G-H).These findings were supported by the bootstrap analysis, based on which, in the whole brain, more voxels were significantly activated when correcting for either vascular covariate compared to applying no correction (Supplementary Figure S3A-B).The two covariates resulted in a similar number of voxels that were significantly activated in response to the task (Supplementary Figure S3C) and in a similar number of voxels in which the association between the task response and the vascular covariate was significant (Supplementary Fig. S3D).However, applying BOLD-CBV as a vascular covariate yielded significantly higher values of the T-statistic in these voxels compared to    applying BOLD-CVR (Supplementary Fig. S4C-D).

Discussion
This study aimed to assess the dependence of the breath-hold BOLD signal on CBV dHb and CVR and whether a BOLD-signal-derived marker of CBV may be a practically useful alternative to BOLD-signal-derived CVR in accounting for between-subject variation in task response modeling in group-level BOLD fMRI analysis.Simulations performed using a Monte Carlo framework highlighted a stronger dependence of the breath-hold-induced BOLD signal change on CBV dHb than the relative change in CBF, that is, the CVR, largely due to the greater variability in CBV dHb .For stimulus-evoked BOLD fMRI, moderate and localized increases in the group-level statistical significance of activation were found when accounting for the individual variability in either BOLD-CBV or BOLD-CVR, with BOLD-CBV resulting in the strongest correction and the largest correlation with the response to the task.

Measuring the deoxyhemoglobin-sensitive CBV
Breath-hold-induced BOLD signal variation in the SSS was similar to the variation of PetCO 2 (Fig. 5) as expected from the vasodilatory effect of hypercapnia in increasing CBF and decreasing OEF.Monte Carlo simulations showed a strong dependence of the breath-hold BOLD fMRI signal change on CBV dHb (Fig. 4) and results in vivo showed a significant correlation across subjects between breath-hold induced BOLD signal oscillations and BOLD-CBV.Thus, the method proposed here for calculating BOLD-CBV based on breath-hold BOLD fMRI offers a practical marker of CBV dHb and highlights that the breath-hold induced BOLD signal variation is, in practice, likely to be more reflective of the greater variation in CBV dHb than that of CVR.
Alternative methods for estimating CBV dHb include venous refocusing for volume estimation (VERVE) (Stefanovic and Pike, 2005), vascular space occupancy (VASO) (Lu et al., 2003), and hyperoxic BOLD fMRI (Bulte et al., 2007).VERVE exploits the dependency of the spin-spin relaxation rate of deoxygenated blood on the refocusing interval.It aims to isolate deoxygenated from oxygenated blood by applying a train of refocusing pulses with variable inter-pulse spacing.In contrast with our proposed method, VERVE provides CBV dHb estimates in a single slice, thus unsuitable for correcting group-level fMRI analyses.Conversely, VASO provides three-dimensional maps of simultaneous CBV and BOLD changes by exploiting differences in T 1 longitudinal relaxation times between blood and tissue.However, VASO provides a global estimate of the arterial plus venous CBV, whereas the stimulus-evoked BOLD response is mainly linked to changes in CBV dHb .Finally, hyperoxic BOLD fMRI measures the change in BOLD signal in response to an O 2 gas challenge, which increases venous oxygen saturation in the mild hyperoxic regime (i.e., with no expected decrease in CBF).It then calculates CBV dHb maps by normalizing the map of the change in CBV relative to a purely venous voxel in the SSS.In the present study, signal normalization was also performed relative to a large vein, albeit no exogenous gas inhalation challenge was required, resulting in the easy implementation of the proposed method.However, it must be noted that our proposed method requires the assumption of uniform fractional variation in CBF (CVR) across the brain during the breath-hold (Taneja et al., 2020;Zhao et al., 2021).

Improving the sensitivity of stimulus-evoked BOLD fMRI
The increased sensitivity of group-level BOLD fMRI when using SRT+BOLD-CBV and SRT+BOLD-CVR agreed with previous results demonstrating that correction for breath-hold responsiveness or vascular reactivity increases the spatial extent and statistical significance of detected activation (Murphy et al., 2011).Furthermore, a positive correlation was observed between the response to the task and each vascular covariate, mainly in the precentral gyrus (Fig. 6), confirming that significant between-subject variability in task-induced BOLD response is associated with individual variation in vascular characteristics.The breath-hold modeling results and the greater group-level correlation between task responses and BOLD-CBV than BOLD-CVR suggest that, in practice, a marker of CBV explains more variation in the BOLD signal response than a BOLD-based marker derived from end-tidal CO 2 measurements does.Thus, the vascular covariate BOLD-CBV captures, better than BOLD-CVR, between-subject variability in group-level BOLD fMRI.
The effects of adding the vascular covariates to the group-level model were subtle but positive in enhancing statistical significance.However, the effects were somewhat heterogeneous between brain regions engaged in the task.We suggest that studies seeking to implement a between-subject correction for vascular characteristics may adopt breath-hold-based maps of BOLD-CBV rather than breath-hold-based maps of BOLD-CVR derived from end-tidal CO 2 measurements.In addition to improving statistical significance, measuring BOLD-CBV requires only estimating the breath-hold-induced signal variation in the SSS.It does not require recordings of end-tidal CO 2 , rendering the procedure much easier to implement experimentally, thus reducing error propagation from potentially inaccurate measurements of endtidal CO 2 to the estimated BOLD-CVR.The reduced susceptibility to experimental errors of breath-hold-based BOLD-CBV compared to breath-hold-based BOLD-CVR may explain why BOLD-CBV most accurately represented the inter-subject physiological variability in response to a task (Figs. 6 and 7), despite BOLD-CVR intrinsically accounting for more sources of physiological variability in the BOLD signal ([Hb], OEF 0 , CBV dHb ) by design (Eq.( 2)).

Limitations
The present study presents some limitations.First, Monte Carlo simulations demonstrated that the BOLD variability following an isometabolic vasodilation is mainly driven by the variability in the deoxyhemoglobin-sensitive CBV.However, once the intrinsic sensitivity to the different physiological parameters is determined by the simulated underlying physical and physiological processes, this result depends on the range of values assumed for the simulation parameters.Although the ranges used in the simulations were plausible in healthy subjects (Chiarelli et al., 2022a;Leenders et al., 1990;Zhao et al., 2021), they may be altered in disease, potentially altering the proportion of BOLD signal variation caused by CBV or CVR variation.Moreover, in Fig. 4E, the slope value was smaller than 1.This result was probably related to the extravascular effect of the BOLD signal affecting the value of BOLD-CBV.Indeed, the extravascular component may account for up to 30% of the BOLD signal regarded as intravascular, thus artificially increasing BOLD-CBV (Boxerman et al., 1995).Further, there was a residual negative dependence of BOLD-CBV on OEF 0 (Fig. 4F) which was possibly caused by the supralinear dependence of the intravascular BOLD signal modulation on oxygen saturation.This nonlinear dependence was induced by two factors.First, the intravascular relaxation rate has a quadratic dependence on oxygen saturation (Eq.7).Second, the intravascular relaxation rate is substantially modulated, making evident its nonlinear contribution to the MRI signal (Eq.8).Future studies may explore how to further model the dependence of the intravascular BOLD signal on oxygen saturation and, hence, reduce the dependence of the normalized BOLD signal on OEF 0 .
A second limitation was due to the BOLD-CBV maps being calculated via normalization of the breath-hold BOLD fMRI images relative to one large vein, assuming low intra-subject spatial variability for CVR.In a healthy subject population, this assumption was supported by observations of CVR having a low spatial variability across the brain (Taneja et al., 2020).In disease, these assumptions would need to be carefully evaluated.
A third technical limitation was linked to the average diameter of the distal SSS, which, in the general adult population, is equal to 4.91 ± 1.20 mm (Durst et al., 2016).Undersampling the SSS ROI through alignment with the breath-hold BOLD fMRI image (resolution = 3.43 × 3.43 × 2.00 mm 3 ) probably induced partial volume effects with the surrounding brain parenchyma.To mitigate this shortcoming, the intravascular signal in the SSS was calculated using an algorithm based on the temporal and spatial characteristics of the breath-hold BOLD fMRI data.Future studies, especially at higher field strengths (e.g., 7 T), could explore employing higher resolutions for the BOLD fMRI acquisitions to mitigate partial volume artifacts and improve breath-hold BOLD signal quantification in the SSS.
Further limitations were linked to the assumptions of isometabolism and constant OEF changes throughout the brain.In BOLD fMRI for physiologic measurement, isometabolism is commonly assumed during mild hypercapnia (Germuska and Wise, 2019).On the one hand, apneic challenges such as breath holding may be slightly pro-metabolic (Rodgers et al., 2013); on the other hand, the ratio between the percentage increase in CMRO 2 and the percentage increase in CBF was observed to be small, within 15% (Rodgers et al., 2013).Hence, it is reasonable to assume isometabolism.By contrast, it may be more appropriate to assume a constant change in OEF or CVR throughout the brain in response to breath-hold in groups of individuals who are unlikely to experience strong regional vascular alterations.These include healthy volunteers (both young and old) as well as patient groups without focal pathologies.
Finally, recommendations on using BOLD-CBV or BOLD-CVR as vascular covariates to account for inter-subject variability in higherlevel stimulus-evoked BOLD fMRI may not hold in patient populations.Indeed, stimulus-evoked BOLD fMRI may exhibit false-negative activation in brain areas with impaired neurovascular coupling (Blicher et al., 2012;Para et al., 2017;Siero et al., 2015;Zacà et al., 2014).Moreover, in the present study, care was taken to realign the temporal dynamics of the BOLD response in the SSS with that of GM tissue, as these compartments may present a relative time shift in the BOLD response.However, in disease, the breath-hold-induced BOLD response of tissue may have a different dynamic from that observed in healthy subjects, induced, for example, by stealing effects resulting from impaired vascular vasodilatory reserve (Sobczyk et al., 2014).Thus, further studies are needed to explore these relationships in specific patient populations.

Conclusions
A greater proportion of variation in blood oxygenation level dependent (BOLD) signal changes arising from breath-hold is explained by variation in cerebral blood volume (CBV) than in cerebrovascular reactivity (CVR).The present study proposes a method for estimating BOLD-CBV as a practical marker for deoxyhemoglobin-sensitive blood volume (CBV dHb ) based on breath-hold BOLD fMRI.This method explains more task-induced BOLD signal variability of vascular origin than BOLD-CVR, is practical to calculate, as it does not require the measurement of end-tidal CO 2 , and is reflective of blood volume, an underlying physiological characteristic of brain vasculature.

Funding
The study was supported by the MS Society UK.

Fig. 2 .
Fig. 2. BOLD-CBV and BOLD-CVR maps.Maps of BOLD-CBV and BOLD-CVR are shown for one representative subject (A-B) and averaged across subjects in MNI space (C-D).The same sagittal slice is displayed on each row.A sagittal view was chosen to highlight the vascular component on these maps.

Fig. 5
Fig.5shows the mean and standard deviation (SD) of the PetCO 2 traces and SSS BH-BOLD signals across subjects during the repeated breathholding.These signals had a similar time course and were strongly correlated (r = 0.84, p < 0.001), indicating that the SSS BH-BOLD signal approximates the form of the end-tidal measurements of exhaled CO 2 , although being expressed in different units.For data acquired in vivo, in the GM the average BOLD-CBV was positively and significantly correlated with the average BOLD signal changes induced by the breath-hold task (r = 0.42, p = 0.03).This result was in line with the simulations shown in Fig.4A.Fig.6shows the group response to the SRT task estimated by the three models (A-L) and the association between the task response and each vascular covariate (M-T).Fig.7reports the group-level statistical distribution of Z scores corresponding to the association between the task response and each vascular covariate in each ROI.The statistical distribution of Z scores corresponding to brain activation in response to

Fig. 4 .
Fig. 4. Monte Carlo simulation results.Scatter plots showing the BOLD or the BOLD-CBV signal (i.e., the BOLD signal normalized to its intravascular component) plotted against each of these variables: the volume of deoxygenated blood (A, E), the baseline oxygen extraction fraction (B, F), the hemoglobin level (C, G) and the fractional change in cerebral blood flow (D, H).A least-squares line is shown superimposed on each scatter plot.Each plot's title also reports the corresponding correlation (r) and p-values.Abbreviations.BOLD: blood oxygenation level-dependent signal; CBF: cerebral blood flow; CBF 0 : baseline CBF; CBV dHb : blood volume responsible for generating the BOLD signal; CVR: cerebrovascular reactivity; Hb: hemoglobin; OEF 0 : baseline oxygen extraction fraction.

Fig. 5 .
Fig. 5. Variability of PetCO 2 and intra-SSS breath-hold BOLD signals.Temporal means and SDs across subjects are shown for the PetCO 2 signal (A) and the breathhold-BOLD signal in the SSS (B).The PetCO 2 trace is shown in mmHg, whereas the SSS trace is shown as the percentage change of the BOLD signal.Abbreviations.BH-BOLD: breath-hold blood oxygenation level-dependent signal; SD: standard deviation; PetCO 2 : partial pressure of end-tidal carbon dioxide; SSS: superior sagittal sinus.

Fig. 6 .
Fig. 6.Group-level fMRI analysis.Z score maps showing the response to the SRT task according to the three models (A-L) and the association between the response to the SRT task and each vascular covariate (M-T) are overlaid on the same representative slices of the MNI template (slices = 40, 65, 70, and 72).Each column shows a different slice.Abbreviations.A/P: anterior/posterior; MNI: Montréal Neurological Institute; R/L: right/left; SRT: serial reaction time.

Fig. 7 .
Fig. 7. Association between task response and vascular covariates.For each ROI, the histograms show the distribution of Z score values corresponding to the association between the response to the SRT task and the vascular covariate in 2-mm isotropic MNI space.A nonparametric kernel density function fitted to the data for each histogram is shown as a continuous line.In each plot and for each distribution, the legend shows the number of voxels in the histogram and the median value of the distribution.Abbreviations.L/R: left/right; mdn: median; SRT: serial reaction time; vox: voxels.
This work was partially conducted under the framework of the Departments of Excellence 2018-2022 initiative of the Italian Ministry of Education, University and Research for the Department of Neuroscience, Imaging and Clinical Sciences (DNISC) of the University of Chieti-Pescara, Italy.This project has received funding from the European Union's Horizon Europe research and innovation programme under the Marie Skłodowska-Curie grant agreement No 101066055 -acronym HERMES.Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Executive Agency (REA).Neither the European Union nor the granting authority can be held responsible for them.Funded in part by the European Union -NextGenerationEU under the National Recovery and Resilience Plan (NRRP), Mission 4 Component 2 -M4C2, Investment 1.5 -Call for tender No. 3277 of 30.12.2021 Italian Ministry of Universities Award Number: ECS0000004, Project Title: "Innovation, digitalisation and sustainability for the diffused economy in Central Italy," Concession Degree No. 1057 of 23.06.2022 adopted by the Italian Ministry of Universities, CUP: D73C22000840006.This project was partially supported by the UK Engineering and Physical Sciences Research Council (EP/S025901/1).This research was funded in part, by the Wellcome Trust [RGW: 104943/Z/14/Z; MG: 220575/Z/20/Z; KM: 200804/Z/16/Z, 224267/ Z/21/Z].For the purpose of open access, the author has applied a CC BY E. Biondetti et al.

Table 2
MRI protocols and acquisition parameters.