Problems and solutions in quantifying cerebrovascular reactivity using BOLD-MRI Open Access

The cerebral vasculature is finely regulated in the healthy brain, such that adequate blood flow is in continuous supply to the cerebral tissue. As such, in addition to responding to changes in perfusion pressure via autoregulation (Lassen, 1959) and neuronal activity via neurovascular coupling (Attwell et al., 2010), the cerebral vasculature is remarkably responsive to CO₂, and to an extent, hypoxia (Kety & Schmidt, 1948; Mardimae et al., 2012). Specifically, it is well known that increased partial pressure of CO₂ (PCO₂) results in transient vasodilation and increased cerebral blood flow (Kety & Schmidt, 1948). Researchers and clinicians have long exploited this phenomenon by administering vascular challenges (e.g., increased CO₂ in the inspired air) while recording the magnitude of the corresponding vascular response (Fisher & Mikulis, 2021; Sleight et al., 2021)—this practice is termed cerebrovascular reactivity (CVR) imaging.

CVR imaging is commonly used to assess the vasodilatory capacity of stenotic vessels and additional downstream vessels in steno-occlusive disease (Duffin et al., 2018; Hartkamp et al., 2017; Sleight et al., 2021; Venkatraghavan et al., 2018). CVR imaging has also been used in various other diseases, including glioma, dementia, and stroke (Fierstra et al., 2018; Krainik et al., 2005; Pillai & Zacá, 2011; Yezhuvath et al., 2012). Given the primary use case of CVR imaging, the magnitude of the vasodilatory response is most directly assessed by measuring the cerebral blood volume (CBV) percent change per partial pressure change in CO₂ (i.e., $ΔPCO2$⁠) in the vessel(s) of interest. This index (⁠$CVRCBV$⁠) can thus be termed the ground-truth CVR^* (even though other definitions are also used in the literature, as will be described):

Although it is the ideal measure, $CVRCBV$ is seldom reported in the literature (Donahue et al., 2009; Lu et al., 2003) due to challenges of imaging changes in CBV. Thus, CVR is more often and easily estimated in terms of the resulting percent change in cerebral blood flow (CBF), and to a lesser extent, blood velocity (v); these physiological indices are hereafter termed $CVRCBF$ and $CVRv$⁠, respectively:

Magnetic resonance imaging (MRI) paired with a vascular challenge is commonly used for CVR imaging and quantification in the human brain (see references in Fisher et al., 2018; Sleight et al., 2021). MRI-quantified CVR represents an average response across many dilating blood vessels. Various MRI-based techniques, not limited to arterial spin labeling (ASL), vascular space occupancy (VASO), blood oxygenation level-dependent (BOLD) MRI, and phase contrast (PC) angiography, have been used for CVR quantification by exploiting changes in either blood volume, flow, or velocity (figure 2 in Sleight et al., 2021). For example, VASO is used to estimate percent change in CBV, ASL is used to estimate percent change in CBF, and PC angiography is used to estimate percent change in v and/or CBF in the large vessels (Donahue et al., 2009; Lu et al., 2003; Taneja et al., 2020; Zhao et al., 2021).

Each CVR method is partially indirect given the confounds of various assumptions and acquisition/analysis parameters. However, BOLD-MRI CVR (CVR_BOLD), which is used to measure T₂^*- or T₂-weighted MRI signal changes associated with CVR-induced blood oxygenation change, arguably provides the most indirect estimation of CVR. In addition, given the relative ease of measurement, CVR_BOLD is implemented in the vast majority—roughly 77%—of CVR-MRI studies (Sleight et al., 2021). To add further complexity, CVR imaging is performed using an array of CO₂ stimulus protocols, which either require the manipulation of $ΔPCO2$ or the harnessing of endogenous $ΔPCO2$ fluctuations. While earlier studies focused predominantly on the former approach (e.g., breath-holding, rebreathing, and prospective end-tidal targeting), resting-state signal-based CVR_BOLD is now becoming increasingly prevalent in the CVR literature (Jahanian et al., 2014; Jahanian et al., 2017; Kannurpatti et al., 2014; Liu et al., 2021; Pinto et al., 2021; Wang et al., 2016; Zang et al., 2007; Zou et al., 2008).

In this paper, we first sought to briefly describe the relationship between the different formulations of physiological CVR (as described in Eqs. 1.1–1.3). We then developed and applied a novel CVR_BOLD simulation framework and interrogated the CVR_BOLD literature to address whether and to what extent CVR_BOLD accurately reflects ground-truth CVR. First, we investigated the relationship between quantified CVR_BOLD and physiological/acquisition parameters, such as the baseline blood oxygenation and magnetic field strength. Using simulations and our own experimental data, we then investigated whether different stimulus protocols, such as breath-holding and resting-state, provide comparable estimates of CVR_BOLD relative to one another. Finally, we simulated example cases of vascular pathology to assess the degree to which CVR_BOLD reflects ground-truth CVR in disease.

As previously described, the ground-truth (i.e., most direct) definition of CVR is $CVRCBV$⁠. However, in the MRI literature, CVR has been described on separate occasions as a CBV change (i.e., in VASO; Donahue et al., 2009; Lu et al., 2003), a v change (i.e., in PC angiography; Hartkamp et al., 2012; Miller et al., 2019), and, most commonly, a CBF change (Cohen & Wang, 2019; De Vis et al., 2015; Hoogeveen et al., 2024; Tancredi et al., 2012; Taneja et al., 2020; Zhou et al., 2015). In addition, CVR has been described as representing vasodilatory capacity in either the large vessels (i.e., using PC angiography; Hartkamp et al., 2012; Miller et al., 2019) or microvasculature. The appropriate definition of CVR then depends on the physiological change of interest, but one must be aware that these measures can provide different representations of vascular physiology and may consequentially provide seemingly contradictory results (see Fig. 1). Thus, great care must be taken when comparing CVR results from imaging methods that exploit different representations of CVR.

Non-invasive CBV-based methods are not available in most sites and are rarely used in the literature for CVR measurements (Lu & Van Zijil, 2012; Sleight et al., 2021); v-based methods do not currently assess microvascular flow due to low spatial resolution in MRI studies (Hartkamp et al., 2012; Miller et al., 2019; Sur et al., 2020; Taneja et al., 2020); CBF-based methods do not usually measure venule/venous signal (due to T₁ decay and water exchange in the capillaries; Le et al., 2012); and may be confounded by the chosen post-labeling delay and blood delivery/transit times (Inoue et al., 2014; Xu et al., 2024). Therefore, although various MRI methods are available for estimating CVR, they are each confounded by limitations in acquisition and availability, such as differences in sensitivity to various blood vessel types. Nevertheless, it is also important to understand how the physiological measures of interest differ from one another in the absence of acquisition-dependent limitations; thus, we first sought to employ simulations to better understand the differences between these ideal physiological measures.

To perform the simulations, the CVR relationships described in Eqs. 1.1–1.3 are reformulated as a function of the change in partial pressure of expired CO₂ (⁠$ΔPETCO2$ (mmHg)), the corresponding global CBF percent change $(ΔCBFp=100·(CBF−CBF0CBF0))$⁠, and a CVR heterogeneity factor (⁠$CVRFactor$⁠). The $CVRFactor$ is a novel term used in this work to represent spatial heterogeneity in physiological CVR across the cerebral vasculature. The $CVRFactor$ represents the factor increase (> 1) or decrease (< 1) in percent flow change for a given simulated voxel relative to the global $ΔCBFp$ in the brain (i.e., if $CVRFactor$ = 2, simulated voxel $CVRCBF$ corresponds to double the $ΔCBFp$⁠, normalized to $ΔPETCO2$⁠). Thus, $CVRFactor$ is effectively a proxy to voxel-specific $CVRCBF$⁠, scaled by $ΔCBFp$ and $ΔPETCO2$⁠.

To quantify $CVRCBV$ (Eq. 2.2) in terms of $ΔCBFp$⁠, Grubb’s relationship is used (Grubb et al., 1974). $CVRCBV$ then depends on $ΔCBFp$⁠, the $CVRFactor$⁠, and Grubb’s exponent (⁠$α$⁠) which governs the relationship between CBF and CBV changes as $(CBVCBV0)=(CBFCBF0)α$⁠.

$α$ was found in earlier studies to be 0.38 (Grubb et al., 1974), but more recent studies suggest it to be roughly 0.29 on average in the human brain (Ito et al., 2003) and 0.18 for the subset of vessels with deoxygenated blood (Chen & Pike, 2010). That is, for the same $ΔCBFp$⁠, veins are expected to yield a smaller ΔCBV relative to capillaries/arteries. Note that this value is also likely to depend on the stimulus/challenge duration (see Uludağ & Blinder 2018).

To quantify $CVRv$ (Eq. 2.3) in terms of $ΔCBFp$⁠, Eqs. 1.3 and 2.1–2.2 were combined with the well-known CBF = v*CBV relationship (here, CBV is %, v is min^-1, and CBF is %/min)^†.

In Figure 1, vessels were simulated as a function of either the $CVRFactor$ or α. For Figure 1A, the simulated vessels for each $CVRFactor$ can be thought of as being from different vascular territories (e.g., middle vs. posterior cerebral arteries); for Figure 1B, the simulated vessels for each $α$ can be thought of as being within the same vascular path (for example, artery to capillary to vein). Given that the literature suggests that $α$ ranges, at the very least, from ~0.18 to 0.4 (Chen & Pike, 2010; Grubb et al., 1974; Ito et al., 2003), we have simulated $α$ just beyond this range in Figure 1B.

As shown in Figure 1A, the chosen definition of CVR is highly determinative of the observed response magnitude; as CBF is proportional to the product of CBV and v, the CBF response should be larger than its constituent parts, which is supported by the literature where average $CVRCBF$ in the human brain is ~3.5–7%/mmHg (Cohen & Wang, 2019; De Vis et al., 2015; Hoogeveen et al., 2024; Tancredi et al., 2012; Taneja et al., 2020; Zhou et al., 2015) and average $CVRCBV$ is ~0.4%/mmHg (when normalized to the expected $ΔPETCO2$ during breath-holding) (Donahue et al., 2009; Lu et al., 2003). While $CVRv$ is similarly expected to be lower than $CVRCBF$⁠, CVR is reported in the large vessels when using PC angiography (~4–8%/mmHg (Miller et al., 2019; Sur et al., 2020; Taneja et al., 2020) in healthy controls) and is, therefore, not necessarily expected to be lower than the ASL-reported $CVRCBF$ values from the microvasculature. Despite differences in the absolute CVR values, relative maps (i.e., relative distribution of CVR values) are expected to agree when $CVRFactor$ is varied, given the linearity in Figure 1A.

The relative distribution of CVR values is dependent on the chosen physiological formulation of CVR when there are variations in $α$ (Fig. 1B). It is known in the literature that $α$ is roughly 0.29 for whole brain circulation (Ito et al., 2003) but 0.18 for the subset of vessels with deoxygenated blood (Chen & Pike, 2010). Thus, it can be assumed that $α$ decreases from arterial circulation to venous circulation, in keeping with the fact that arterial diameter is controlled by smooth muscle cells while venous diameter is more controlled by the intravascular pressure change generated at the arterial/capillary level. Logically, $CVRCBV$⁠, $CVRCBF$⁠, and $CVRv$ will be discordant when comparing vessels with different $α$ values (i.e., a voxel with a higher $α$ value (arterially dominated) will yield a higher $CVRCBV$ and lower $CVRv$ than in a voxel with a lower $α$ value (venous dominated), even though $ΔCBFp$ may be identical across these vessels, if for example, they are all within the same vascular path). Therefore, $CVRCBF$ and $CVRv$ will likely not adequately represent regional variations in the ground-truth (i.e., $CVRCBV$⁠) when $α$ is varied.

In summary, the various physiological formulations of CVR cannot be expected to yield the same absolute values (Fig. 1A). In addition, the formulations of CVR should not yield the same relative maps in the presence of regional variations in the types of blood vessels in a voxel (Fig. 1B).

This section was dedicated to understanding how various formulations of physiological CVR (independent of MRI-associated acquisition confounds) differ from one another in both absolute and relative terms. For the remainder of this work, CVR_BOLD, an index most widely utilized for CVR measurements, albeit beholden to the confounds of BOLD-MRI, will be the focus of our discussion.

BOLD-MRI is used to measure T₂- or T₂*-weighted signal changes resulting from regional changes in magnetic field inhomogeneity, predominantly assumed to arise from changes in paramagnetic deoxyhemoglobin (dOHb) (Bandettini et al., 1992; Kwong et al., 1992; Ogawa et al., 1990, 1992) (note that signal changes can additionally arise from concurrent blood volume changes; Bright et al., 2014; Schulman et al., 2024; Thomas et al., 2013; Uludağ, 2010).

CVR_BOLD is typically measured as a percent change in T₂*-weighted signal normalized to $ΔPETCO2$ (Fisher et al., 2018; Liu et al., 2019):

Here, $S0$ is the baseline BOLD signal (prior to the CO₂ challenge) and $SMax$ is the maximum BOLD signal in response to the CO₂ challenge. CVR_BOLD is assumed to reflect CVR, particularly $CVRCBF$⁠, given that increased CBF leads to a decrease in deoxyhemoglobin (dOHb) in the capillaries and veins, and thus, an increase in T₂*-weighted MRI signal (as summarized in Liu et al., 2019 and Uludağ et al., 2009). However, to link $CVRCBF$ with CVR_BOLD, the contribution of various physiological and physical variables must be considered and assessed.

Using a simulation framework, we sought to determine the influence that various parameters (e.g., echo time, field strength, blood volume) have on CVR_BOLD quantification. The simulation framework in this work is adapted from a DSC signal model (Schulman et al., 2022, 2023) informed by prior fMRI/DSC models (Kjolby et al., 2006, 2009; Uludağ et al., 2009).

T₂*-weighted MRI signal (⁠$S$⁠) is a summation of extravascular (EES) and intravascular (IVS) signal contributions, $SEES$ and $SIVS$⁠, respectively (Uludağ et al., 2009):

The index i denotes the specific vascular component within the voxel being simulated (capillary, venule, or arteriole; see Section 2.2.2). Each signal contribution is scaled by the total baseline blood volume relative to the tissue volume (⁠$CBV0$ (%)), each vessel’s volume fraction (⁠$Vi$, reflecting the fraction of blood volume attributable to a given vessel type), and each vessel’s hypercapnia-induced decimal change in blood volume (⁠$ΔCBVi$) (see Section 2.2.4 for details). Note that we assume that the spin density factor (ϕ; the ratio of IVS-to-EES spin density) is equal to 1 (Lu et al., 2003; Uludağ et al., 2009). Therefore, $S$ in our simulations does not consider spin density effects, but this can easily be incorporated in the future by modifying ϕ. In addition, $S$ is normalized such that at an echo time of 0 ms, $S$ is equal to 1.

The signal components are separable as follows:

$R2,con,EES,i*$ and $R2,con,IVS*$ (s^-1) are the extravascular and intravascular magnetic relaxation rates induced by a susceptibility agent, respectively. $R2,0,EES*$ and $R2,0,IVS*$ (s^-1) are the baseline magnetic relaxation rates without a susceptibility agent (i.e., when hemoglobin is fully saturated with oxygen). TE is the echo time (ms), and although it is simulated across a range of values in Section 3.1, TE is set to 30 ms for most of the simulations.

Table 1 summarizes the extravascular and intravascular relaxation rates across field strength found in the literature (Uludağ et al., 2009; Zhao et al., 2007).

The baseline voxel parameters used in the simulations are based on previously reported physiological estimates (Lu & Ge, 2008; Uludağ et al, 2009; Vovenko, 1999) and are summarized in Table 2. However, note that in Sections 3 and 4, these parameters are varied across a range of values to determine their effect on CVR_BOLD quantification.

Tissue blood oxygenation increases during hypercapnia (Jain et al., 2011; Kety & Schmidt, 1948; Poublanc et al., 2021). The change in oxygen saturation resulting from hypercapnia (⁠$ΔYhc$) can be estimated from the well-known cerebral metabolic rate of oxygen (CMRO₂) and CBF relationship (Mintun et al., 1984)—note that CMRO₂ remains relatively constant during mild/moderate hypercapnic challenges (Blockley et al., 2013; Chen & Pike, 2010; Jain et al., 2011; Vestergaard & Larsson, 2019; Zappe et al., 2008), allowing for the derivation of Eq. 6.1. $ΔYhc$ depends on the simulated vessel. $ΔYhc$ for vein/venule (⁠$ΔYhc,v$⁠) depends on $ΔCBFp$⁠, the oxygenation extraction (baseline arteriole oxygenation (⁠$Y0,a$⁠)—baseline venule oxygenation (⁠$Y0,v$⁠)), and the $CVRFactor$ in simulated tissue. To calculate $ΔCBFp$⁠, we first set $ΔPETCO2$ to 9.25 mmHg based on typical maximal changes during $ΔPETCO2$-based CVR studies (Bright & Murphy, 2013; Poublanc et al., 2015; Sleight et al., 2021). We then used literature data/equations (Grüne et al., 2015; Poulin et al., 1996; Ramsay et al., 1993; Sato et al., 2012; Tancredi & Hoge, 2013) to calculate $ΔCBFp$ in response to a $ΔPETCO2$ of 9.25 mmHg (⁠$ΔCBFp$ of ~40%). As previously discussed, $CVRFactor$ is a novel term used in this work to represent spatial heterogeneity in $CVRCBF$ across the cerebral vasculature (i.e., simulated voxels).

$ΔYhc$ for artery/arteriole (⁠$ΔYhc,a$⁠) is set directly to 0, given that in healthy subjects there is negligible oxygen exchange in the arterial/arteriolar segment of the vascular network (Jain et al., 2011; Vu et al., 2023).

$ΔYhc$ for capillary (⁠$ΔYhc,c$⁠) is half of that in simulated venule (based on the fact that the oxygen extraction in capillary (relative to arteriole) is assumed to be approximately half of that in venule (relative to arteriole)).

The total change in oxygen saturation (⁠$ΔY$⁠) for a given vessel can then be described as the summation of its $ΔYhc$ and the oxygenation change resulting from hypoxia (⁠$ΔYhyp$⁠).

$ΔYhyp$ is assumed to be 0 for all simulations, except for breath-holding (see Section 5).

As previously mentioned (Eq. 4), hypercapnia-induced vasodilation (⁠$ΔCBV$⁠) must be included in the CVR modeling (Kety & Schmidt, 1948; Liu et al., 2019). $ΔCBV$ depends on the $ΔCBFp$⁠, $CVRFactor$⁠, and Grubb’s exponent (⁠$α$⁠) (Grubb et al., 1974) in simulated tissue:

As previously noted, $α$ for whole brain has been found to be roughly 0.29 (Ito et al., 2003) but 0.18 for the subset of vessels with deoxygenated blood (Chen & Pike, 2010). Thus, for our simulations we set $α$ to 0.18 for venules, 0.29 for capillaries, and extrapolated to 0.4 for arterioles. Note that $ΔCBV$ is relative (e.g., if $ΔCBV$ = 1.08 (i.e., an 8% change) and $CBV0$ = 4% (typical values for gray matter), CBV would increase to 4.32%).

According to a recent systematic review (Sleight et al., 2021), the vast majority of papers report CVR_BOLD as percent signal change normalized to $ΔPETCO2$ (⁠$CVRBOLD,S$⁠), described in Eq. 3, whereas only 0.4% of studies report CVR_BOLD as a relaxation rate change normalized to $ΔPETCO2$ (⁠$CVRBOLD,R$⁠):

As a reminder, $S0$ is the baseline signal (when $ΔCBFp$ is set to 0), $SMax$ is the maximum signal (e.g., when $ΔCBFp$ is set to 40%), and $ΔR2*$ is the corresponding change in relaxation rate. The preponderance of $CVRBOLD,S$ relative to $CVRBOLD,R$ is very surprising given the well-known relationship between echo time and T₂* signal change (Eqs. 5.1 and 5.2). That is, higher TE yields higher calculated values of $CVRBOLD,S$ (e.g., see figure 1 in Havlicek et al., 2017), whereas $CVRBOLD,R$ should be (almost) TE independent. As previously mentioned, $ΔPETCO2$ is set to 9.25 mmHg for these simulations.

In addition to investigating the effect of TE on CVR_BOLD using our simulations, we surveyed gray matter (GM) and whole brain (WB) $CVRBOLD,S$ values reported in the literature across a range of TEs at 3 T (Fig. 2).

Literature values at 3 T (Atwi et al., 2019; Burley et al., 2021; Dengel et al., 2017; Donahue et al., 2014; Hou et al., 2020; Leung, Duffin, et al., 2016; Leung, Kim, et al., 2016; Liu, Welch, et al., 2017; Ravi et al., 2015, 2016; Sleight et al., 2023; Sobczyk et al., 2016, 2021; Triantafyllou et al., 2011; Van Niftrik et al., 2018; Vu et al., 2023; Zhou et al., 2015) suggest that a positive, linear TE-dependency exists when reporting as $CVRBOLD,S$ (Fig. 2A), in agreement with our simulations (Fig. 2B). In fact, Triantafyllou et al. specifically investigated the dependency of $CVRBOLD,S$ on TE and observed a strong positive linear correlation (Triantafyllou et al., 2011). Note that this relationship assumes that signal is above the noise floor (i.e., we would not necessarily expect this linear relationship to hold true in regions where baseline signal is low and dominated by noise). While there are insufficient data for $CVRBOLD,R$ across a range of TEs in the literature, our simulations show that the TE-dependency is largely eliminated when CVR is calculated as $CVRBOLD,R$ (i.e., $CVRBOLD,R$ values differ by only ~2–4% between TE = 20 ms and TE = 40 ms) (Fig. 2C).

Given these insights, we recommend that researchers and clinicians implement Eqs. 9.1–9.2 for CVR_BOLD,R calculation, to reduce the TE-dependency associated with $CVRBOLD,S$⁠, allowing for an easier comparison of CVR_BOLD values across studies.

An increased contrast-to-noise ratio makes a compelling case for the implementation of increased magnetic field strength in CVR_BOLD studies; however, the relationship between field strength and both $CVRBOLD,S$ and $CVRBOLD,R$ must be considered.

Triantafyllou et al. investigated $CVRBOLD,S$ in GM at 1.5 T and 3 T, reporting an approximate 1.76-fold increase in percent signal change normalized to flow change as field strength doubled (Triantafyllou et al., 2011). In a breath-holding experiment, Peng et al. found a 1.66-fold increase in percent signal change as field strength doubled from 1.5 T to 3 T (Peng et al., 2020). Driver et al. investigated $CVRBOLD,R$ in GM at 3 T and 7 T, reporting an approximate doubling in relaxation rate change as field strength increased by 2.33-fold (Driver et al., 2010). Note that in comparison to CVR_BOLD, ASL-derived CVR (CVR_ASL), as expected, does not significantly depend on echo time or field strength (Nöth et al., 2006). As these studies only compared results across two field strengths, it is difficult to identify whether the field strength relationship is linear across a wider range. Therefore, we used our simulations to determine the relationship between field strength (simulated across 1.5 T, 3 T, 4.7 T, and 7 T; see Table 1) and quantified $CVRBOLD,S$ and $CVRBOLD,R$ (Fig. 3).

GM $CVRBOLD,S$ values reported in the literature (Atwi et al., 2019; Bhogal et al., 2014; Bhogal et al., 2015; Blockley et al., 2011; Burley et al., 2021; Dengel et al., 2017; Donahue et al., 2014; Kassner et al., 2010; Leung, Duffin, et al., 2016; Leung, Kim, et al., 2016; Liu, Welch, et al., 2017; Sleight et al., 2023; Sobczyk et al., 2016, 2021; Stefanovic et al., 2006; Vesely et al. 2001; Vu et al., 2023; Zhou et al., 2015) empirically demonstrate a logarithmic relationship (R² = 0.54) with field strength (Fig. 3A), although more data are likely needed to verify this relationship.

Like the aforementioned literature values, the simulations show that CVR_BOLD does not increase linearly, but rather in a quasi-logarithmic manner. In addition to a global scaling effect, simulated voxels with higher CBV yield disproportionately larger increases in $CVRBOLD,R$ and $CVRBOLD,S$ as a function of field strength—in theory, this would likely result in regional scaling differences between CVR_BOLD maps of different field strengths in regions containing large veins/arteries. However, in comparing regions with low blood volume (i.e., majority of the GM and WM), regional scaling differences are minimal.

Given the relationship’s deviation from linearity, a non-linear correction factor can be implemented to remove the influence of field strength on CVR_BOLD measures. Upon fitting the $CVRBOLD,R$ and $CVRBOLD,S$ simulations, we found that logarithmic-radical fits best corrected for $CVRBOLD,R$ (R² = 0.999) and $CVRBOLD,S$ (R² = 0.99):

Therefore, for inter-field strength comparisons of CVR values obtained from similar populations, we recommend that a field-dependent correction factor, like Eqs. 10.1–10.2, is utilized prior to reporting CVR_BOLD values to ensure better inter-field strength comparability. Ideally, this correction factor would be derived and validated empirically by collecting CVR_BOLD data on the same subjects at multiple magnetic field strengths and then compared with the relationship in the simulations.

Note that for simplicity, we only quantifyCVR_BOLDas $CVRBOLD,R$ for the remainder of this work. Although not explicitly shown, $CVRBOLD,S$ results are qualitatively identical to $CVRBOLD,R$ in the forthcoming relationships.

The quantitative relationships between CVR_BOLD and various tissue/vascular properties, such as the hematocrit (Hct) and arteriolar blood volume, remain understudied. We explored the literature and used our simulations to understand the dependency of CVR_BOLD on a comprehensive set of physiological parameters.

One parameter that is expected to influence CVR_BOLD is the underlying baseline CBV (CBV₀), in agreement with the literature (Blockley et al., 2013; Buxton et al., 2004; Davis et al., 1998; Yablonskiy & Haacke 1994), where the relationship between CBV₀ and the percent BOLD change (and/or $ΔR2*$⁠) has been shown to be linear (specifically, see Eq. 1 in Davis et al., 1998 and Eq. 9 in Buxton et al., 2004); this is because the change in the absolute amount of dOHb (susceptibility agent) in a voxel during hypercapnia is scaled by CBV₀, which ultimately scales the observed BOLD signal change. Thus, according to this relationship derived from these standard BOLD signal models (Buxton et al., 2004; Davis et al., 1998), doubling CBV₀ results in a doubling of the measured BOLD-MRI signal change, and thus, CVR_BOLD. Consequently, comparing CVR_BOLD,GM with CVR_BOLD,WM, for example, is heavily biased by the roughly doubled CBV₀ in GM (CBV₀ ~ 4%) relative to WM (CBV₀ ~ 2%). Therefore, making the claim that “CVR_BOLD reflects $CVRCBF$” is just as valid as claiming that “CVR_BOLD reflects CBV_0,” which poses a considerable limitation on CVR_BOLD in its current state. To clarify, if CVR_BOLD reflected CVR_CBV (that is, the change in CBV), this would be a benefit, not a limitation. However, CVR_BOLD scales linearly with baselineCBV(CBV₀). This is a limitation as it confounds the reactivity measure of interest with a baseline measure that is not reflective of the percent CBV change in response to a vascular stimulus.

In agreement with the theory (Blockley et al., 2013; Buxton et al., 2004; Davis et al., 1998; Yablonskiy & Haacke 1994), CVR_BOLD,GM has empirically been found to be roughly double CVR_BOLD,WM in the literature, while CVR_ASL,GM has been found to be roughly the same as CVR_ASL,WM (Taneja et al., 2020; Xu et al., 2024). Since DSC-MRI calculated CBV₀ (Ostergaard, 2005; Schulman et al., 2023) and CVR_BOLD both scale roughly linearly with CBV₀, dividing CVR_BOLD by CBV₀ (for example, measured separately using DSC) is a promising solution for removing the confound of CBV₀ from measures of CVR_BOLD. The DSC data can be acquired using Gd-based contrast (in clinical scenarios, this is often routinely administered for glioma and neurovascular disease) or dOHb-based contrast (administered using sequential gas delivery systems, as is used in many CVR studies to induce hypercapnia). This novel approach could then be validated by comparing with CVR_ASL which is assumed to be independent of CBV₀ (Petersen et al., 2006; Taneja et al., 2020).

The observed tissue relaxation rate change in CVR_BOLD is a result of increased blood flow, reduced oxygen extraction, and ultimately, a lower proportion of deoxygenated hemoglobin in the capillaries and veins. In arterial/arteriolar blood where oxygenation is near full saturation (Severinghaus, 1979) and where oxygen does not generally exchange with tissue, oxygenation-induced relaxation rate change is negligeable. We used our simulations to demonstrate the effect of the simulated proportion of arteriolar CBV₀ (CBV_A (%)) on the measured CVR_BOLD (Fig. 4A).

While $CVRCBF$ and $CVRCBV$ are held constant, CVR_BOLD reduces drastically as the simulated tissue’s vessels become arterially dominated (theoretically, there is no transverse relaxation rate change when the voxel only contains arterial/arteriolar blood). Therefore, voxels with a low CVR_BOLD may have high vasodilatory capacity and only appear to have a low CVR due to a higher proportion of fully oxygenated blood. In other words, CVR_BOLD has a low sensitivity to CVR in arteries/arterioles (see Schulman et al., 2024, which shows how CVR_BOLD in arteries can even be negative due to the displacement of cerebrospinal fluid). There is no straightforward way to mitigate this shortcoming, but studies should acknowledge the lack of oxygenation-induced arterial/arteriolar relaxation change as a limitation of CVR_BOLD.

While baseline arterial/arteriolar blood oxygenation (⁠$Y0,A$ (%)) is usually ≥ 97% in humans (Chan et al., 2013), it can be markedly lower (within the range of 80–90%) in certain patient populations, such as those with pulmonary disease (Echevarria et al., 2021). Thus, we sought to investigate the influence of $Y0,A$ on CVR_BOLD (Fig. 4B). Of note, it is known that reducing $Y0,A$ beyond roughly 80% (i.e., ~50 mmHg) results in a hypoxia-driven cerebrovascular response (Mardimae et al. 2012); therefore, we constrained the $Y0,A$ parameter to avoid simulating this additional hypoxia-mediated increase in CBF.

The relationship between $Y0,A$ and CVR_BOLD is heavily dependent on the CBV₀ of the simulated voxel (i.e., CVR_BOLD magnitude increases non-linearly with $Y0,A$ in tissue (CBV₀ ~2–8%) but decreases non-linearly with $Y0,A$ in highly vascular voxels (e.g., CBV₀ of 32%)). The observation for the highly vascular voxels is in agreement with our previous DSC MRI work, which found that $∫ΔBOLD(t)$ induced by susceptibility agent was negatively correlated with $Y0,A$ in voxels with larger vessels (Schulman et al., 2023), owed to the quadratic dependency of the intravascular relaxation rate on blood oxygenation. However, the observed CVR_BOLD increase as a function of $Y0,A$ in tissue is owed to a larger overall extravascular relaxation rate change magnitude at higher $Y0,A$ as a result of concurrent ∆CBV-mediated relaxation rate change (opposite in sign to the CVR-mediated relaxation rate change), which is larger in magnitude at lower $Y0,A$ (refer to Schulman et al., 2024). In general, these simulations reveal a complex relationship between the measured CVR_BOLD and $Y0,A$ which is important for consideration in hypoxemic patients. In future work, incorporating the measured $Y0,A$ (estimated from a patient’s S_pO₂) into tissue CVR_BOLD modeling may bypass this dependency in patients with abnormally low $Y0,A$ values.

While Hct is measured to be around 40–54% for men and 36–48% for women (Billett, 1990), Hct can be as low as ~20% in anemic patients (Afzali-Hashemi et al., 2021; Kuwabara et al., 2002; Nur et al., 2009). While the literature suggests that certain forms of anemia may result in reduced CVR_BOLD (Kosinski et al., 2017), it is unclear to what extent Hct attenuates CVR_BOLD when $CVRCBF$ and $CVRCBV$ are unchanged. To examine this relationship, we simulated (see Section 2.2.1 and Table 1) CVR_BOLD for control (Hct = 44%) and anemic (Hct = 21%) tissue/blood (Fig. 4C). For these simulations, we assumed that reduced Hct yields a disproportionally large increase in baseline blood flow, resulting in halving of the baseline OEF (as observed in Vu et al., 2017). Although not shown, simulating a doubling of the baseline OEF or maintenance of the baseline OEF yields the same pattern observed in Figure 4C relative to control, albeit with different amplitudes.

According to the simulations (Fig. 4C), even in the absence of any changes to the underlying CVR (i.e., $CVRFactor$ = 1), the quantified CVR_BOLD is reduced in simulated anemia by ~40%, relative to normal tissue—this appears regardless of whether the $CVRFactor$ is low (0.5) or high (2). These results have important implications when assessing CVR_BOLD in anemic patients, and even when comparing results across sex—namely, while reduced $CVRCBF$ might partially account for the observed CVR_BOLD reduction in anemia relative to controls (Afzali-Hashemi et al., 2021; Kosinski et al., 2017; Kuwabara et al., 2002; Nur et al., 2009; Sayin et al., 2022) and in women relative to men (Kassner et al., 2010), the simulations indicate that CVR_BOLD increases as a function of Hct even when physiological CVR (i.e., $CVRCBV$ or $CVRCBF$⁠) is unchanged. Upon validating this dependency experimentally (i.e., observing CVR_BOLD as a function of Hct in healthy subjects) and simulating across a larger range of Hct values, incorporating an individual’s Hct into a Hct-based correction factor may allow for the correction of this CVR_BOLD dependency.

Note that variations in baseline CMRO₂ would theoretically confound the relationship between Hct and CVR_BOLD as it would alter $Y0$ (i.e., increased baseline CMRO₂ would yield a decrease in $Y0$⁠). CVR_BOLD is dependent on $Y0$⁠, as we have shown above (Fig. 4B).

Breath-holding is a well-accepted method in CVR_BOLD studies (Bright & Murphy, 2013; Dlamini et al., 2018; Fierstra et al., 2013; Peng et al., 2020; Ratnatunga & Adiseshiah, 1990); however, there are notable limitations. One of the major limitations of breath-holding in comparison with controlled hypercapnia methods (i.e., dynamic end-tidal forcing or prospective end-tidal targeting) is the confounding effect of concurrent hypoxia (Tancredi & Hoge, 2013) whose level is dependent on the breath-hold duration (Sasse et al., 1996). To date, it is unclear how the degree of hypoxia propagates error into the quantified CVR_BOLD values.

Using hypoxia and hypercapnia values reported in previous works (refer to Supplementary Fig. 1 and Supplementary Table 1), we simulated CVR_BOLD as a function of the breath-hold duration. To reiterate, breath-holding (Fig. 5B) is associated with hypoxia whereas controlled hypercapnia (Fig. 5A) is not (here, breath-hold duration is just a proxy to the $ΔPETCO2$ stimulus magnitude; see Supplementary Table 1). In addition, while some researchers (Bright & Murphy, 2013; Golestani et al., 2016; Lipp et al., 2015) normalize breath-hold-acquired CVR_BOLD values to $ΔPETCO2$ (red), others do not (blue) (Cohen & Wang, 2019; Jahanian et al., 2017). Thus, we investigated the quantitative consequences of these different analysis and hypercapnic stimulus design choices.

In the controlled hypercapnia simulations, $ΔPETCO2$-normalized CVR_BOLD progressively decreases, though minimally, and $ΔPETCO2$-unnormalized CVR_BOLD progressively increases with the stimulus magnitude until $ΔPETCO2$ saturates after changing by ~10.5 mmHg (Sasse et al., 1996) (Fig. 5A). The mild decrease in normalized CVR_BOLD is initially perplexing, especially given that $CVRCBF$ mildly increases as a function of $ΔPETCO2$ (Supplementary Table 2). This becomes clearer when recognizing that BOLD begins to saturate with progressively larger increases in $ΔPETCO2$ and $ΔCBFp$⁠, as has been observed in the literature (Bhogal et al., 2014, 2015; Duffin et al., 2017; Hoge et al., 1999) and in agreement with Eq. 6.1.

In the breath-holding simulations, $ΔPETCO2$-normalized CVR_BOLD decreases substantially, and $ΔPETCO2$-unnormalized CVR_BOLD increases with breath-hold durations up to 20 s and then decreases as breath-hold duration increases further (Fig. 5B). This counterintuitive decline is due to a competing hypoxic effect (i.e., increased paramagnetic dOHb during hypoxia results in a positive relaxation rate change (Poublanc et al., 2021; Sayin et al., 2023; Schulman et al., 2023; Vu et al., 2021), which counteracts the negative relaxation rate change from CVR); the summation is then either an increase or decrease in relaxation, depending on the larger contributor (see Supplementary Fig. 2). The finding that $ΔPETCO2$-unnormalized CVR_BOLD values are more consistent at moderate breath-hold durations is observed in previous work, where 10, 15, and 20 s breath-hold durations were compared for CVR_BOLD measurement (Bright & Murphy, 2013); here, CVR_BOLD values were found to be relatively similar between 15 s and 20 s (difference of ~0–10%) as opposed to between 10 s and 15 s (difference of ~50–60%), similar to our simulated findings.

Note that studies have reported different $ΔPETCO2$ values associated with breath-hold durations (e.g., Sasse et al., 1996 vs. Bright & Murphy, 2013). For this paper, we opted to implement values from Sasse et al. (1996) as they span a large range of breath-hold durations and were collected from arterial blood. However, if we input the $ΔPETCO2$ values measured from Bright and Murphy (2013) into our model (instead of those from Sasse et al., 1996), we obtain $ΔPETCO2$-normalized CVR_BOLD values that are much more stable across breath-hold duration (like in Bright & Murphy, 2013): ~0.057, 0.055, and 0.051 s^-1/mmHg for 10, 15, and 20 s breath-holds, respectively.

Ultimately, breath-holding is expected to yield an underestimation of CVR_BOLD (normalized and unnormalized), particularly for longer breath-hold durations, in comparison with controlled hypercapnia due to the associated hypoxic effect (Supplementary Table 1). This is in keeping with the literature where CVR_BOLD values in GM tend to be ~16% lower, on average, for breath-holding versus controlled hypercapnia (Atwi et al., 2019; Bright & Murphy, 2013; Dengel et al., 2017; Donahue et al., 2014; Hou et al., 2020; Leung, Duffin, et al., 2016; Leung, Kim, et al., 2016; Lipp et al., 2015; Liu, Welch, et al., 2017; Moia et al., 2020, 2021; Murphy et al., 2011; Pinto et al., 2016; Ravi et al., 2015; Sleight et al., 2023; Sobczyk et al., 2016; Van Niftrik et al., 2018; Zhou et al., 2015). Additionally, based on the simulations, breath-hold durations of ~20–25 s maximize the contrast-to-noise ratio in breath-hold CVR_BOLD experiments. However, the simulation results suggest that those seeking to limit the hypoxic effect (i.e., those seeking breath-hold CVR_BOLD values in agreement with controlled-hypercapnia CVR_BOLD values) should aim for a breath-hold duration of ~10–15 s—this range is anyway expected to be more translatable to patients.

Over the past two decades, researchers have claimed that CVR_BOLD can be measured using low-frequency resting-state oscillations (Jahanian et al., 2014; Kannurpatti et al., 2014; Wang et al., 2016; Zang et al., 2007; Zou et al., 2008), with the rationale being that these fluctuations are temporally correlated with spontaneous $PETCO2$ fluctuations (Wise et al., 2004). In general, resting-state CVR_BOLD can largely be grouped into signal variation methods, where signal fluctuation magnitude is considered, and regression methods, where magnitude and delay are both considered with respect to a regressor; these methods have been described at length in previous works (Golestani et al., 2016; Jahanian et al., 2014, 2017; Lipp et al., 2015; Liu, Li, et al., 2017; Pinto et al., 2021; Poublanc et al., 2015; Zou et al., 2008; Zuo et al., 2010).

Although resting-state CVR_BOLD has been shown to correlate well with controlled hypercapnia and breath-holding methods in a few studies (Jahanian et al., 2017; Liu et al., 2021), its validity as a technique for measuring CVR is unclear. In particular, the measured fluctuation magnitude, even when filtered to maximally correlate with $PETCO2$ fluctuations (Liu, Li, et al., 2017), is likely confounded by other processes which do not reflect CVR—most notably, neuronal activity-induced BOLD fluctuations (Smitha et al., 2017; Van den Heuvel and Hulshoff Pol, 2010).

We conducted a voxel-wise statistical comparison between breath-hold and resting-state methods (using our own unpublished 3 T and 7 T data) to determine whether resting-state is systematically different from breath-hold CVR_BOLD (denoted in the experimental data as CVR*), and if so, where in the brain these differences are most pronounced. The quantification of CVR* will now be briefly described (see Section 1 of the Supplementary Materials for details).

Using a recently developed method to extract an arterial vasodilatory time course in response to hypercapnia (Schulman et al., 2024), we implemented a novel input regressor for both breath-hold and resting-state paradigms. The resting-state data were bandpass filtered from 0.01 to 0.045 Hz in order to isolate the frequency range of predominate resting $PETCO2$ fluctuations (Wise et al., 2004). For both resting-state and breath-hold paradigms, the relaxation rate time courses were linearly interpolated to a temporal resolution of 0.5 s. Time courses within the arterial regressor mask were multiplied by -1 and averaged voxel-wise to generate the arterial input regressor function (⁠$AIF$⁠). This vertical flipping of the AIF is necessary because the AIF measured in response to hypercapnia is opposite in sign relative to the tissue time courses (see Schulman et al., 2024). The $AIF$ was fit to each tissue voxel by first shifting the $AIF$ forward in 0.5 s steps, up to a maximum of 7 s, and the correlation between the $AIF$ and tissue time courses was recorded for each shift. The temporal shift resulting in the highest correlation was recorded as the voxel’s delay (d). The $AIF$ was then fit to the time course in each voxel (⁠$Tissue(t)$⁠), using the voxel’s delay as an input parameter, assuming the following model:

Here, $CVRBOLD,R*$ is the scaling factor that minimizes least squared error between the tissue and $AIF$ relaxation rate time courses. In terms of the resulting maps, $CVRBOLD,R*$ is effectively a proxy to a regression-based estimation of $CVRBOLD,R$ (the absolute value of $CVRBOLD,R*$ will be different from $CVRBOLD,R$ as values are not normalized to $ΔPETCO2$⁠, but comparing absolute values is not the focus of our analysis). Note that $CVRBOLD,R*$ is not a proxy to steady-state$CVRBOLD,R$⁠, which requires the additional modeling of a hemodynamic response function (HRF) (Poublanc et al., 2015). Finally, $CVRBOLD,R*$ was normalized to the average GM $CVRBOLD,R*$ to obtain the relative CVR maps (⁠$CVR*$⁠). $s$ represents the vertical shift needed to account for baseline MRI signal differences and $ε(t)$ represents the residual error.

In Figure 6, the subject-averaged $CVR*$ maps (which are normalized to average GM) show agreement between field strengths in most regions (i.e., patterns of $CVR*$ values across brain regions are the same at 3 T and 7 T). Note that the same subjects were measured at both field strengths. Correlation values between the AIF regressor and voxel-wise data increase significantly from 3 T to 7 T, particularly for the breath-hold paradigm (breath-hold: p < 0.0005; resting-state: p < 0.005). In Supplementary Figure 4, the delay maps show regional agreement between breath-hold and resting-state maps. At both 3 T and 7 T, there are substantial and significant (p< 0.05) voxel-wise differences between breath-hold and resting-state $CVR*$ maps (Fig. 6B). In comparison with breath-holding, resting-state yields higher $CVR*$ in the visual cortex, pre-central gyrus, and post-central gyrus, but lower $CVR*$ in the subcortical GM (e.g., putamen).

Region-specific resting-state neuronal fluctuations in the visual, motor, sensory, and default mode networks, to name a few, which are active during the resting-state (Gohel & Biswal, 2015) and similar in frequency (Buzsáki & Draguhn, 2004; Penttonen & Buzsáki, 2003; Smitha et al., 2017; Van den Heuvel and Hulshoff Pol, 2010) to CO₂-induced relaxation fluctuations (Wise et al., 2004), are the most likely explanation for these differences between breath-hold and resting-state $CVR*$ maps. While these fluctuations are generally averaged out during the analysis of the breath-hold data and are smaller in amplitude than the breath-hold-induced signal changes, they are generally maintained during bandpass filtering (0.01–0.045 Hz), making it difficult to remove these contributions from resting-state $CVR*$⁠. To make matters more complex, resting-state neuronal fluctuations are highly variable between individuals and depend on their cognitive state during the scan (Gonzalez-Castillo et al., 2021; Han et al., 2023). As observed in Supplementary Figure 5, our resting-state $CVR*$ maps differ substantially between participants, each incorporating aspects of various resting-state networks (particularly the visual network, which is potentially due to the variability of subjects with eyes open vs. closed). The observed inter-subject variability relative to the breath-hold technique limits the generalizability of resting-state CVR_BOLD. Although voxel-wise resting-state CVR_BOLD data have not been widely investigated, it can be observed in Figure 4 from Liu et al. that many subjects displayed regional differences in resting-state versus CO₂-based CVR_BOLD maps (i.e., regional correlations below 0.7), with a similar observation of increased resting-state CVR_BOLD in the visual cortex relative to traditional CVR_BOLD (Liu et al., 2021). Although not calculated here, we also expect differences between breath-hold and resting-state $CVR*$absolute values due to factors such as arterial compliance (Prokopiou et al., 2019).

As a final note, breath-holding-based CVR_BOLD may also be confounded by activations in the brain’s respiratory centers and motor regions that are synchronized with the breath holds. In fact, a recent study found that significant BOLD activations were observed in the pontine respiratory group and raphe nuclei during breath-holding (Ciumas et al., 2023). Thus, although to a lesser extent, breath-holding-based CVR_BOLD is similarly confounded by neuronal activity in certain brain regions.

To better understand the influence of concomitant neuronal activity on CVR_BOLD quantification, we performed a very simple simulation (independent of the simulation framework in Section 2.2). Here, a tissue time course was simulated as the summation of an input CO₂ time course (⁠$CO2(t)$⁠, simple sinusoid with a frequency of 0.033 Hz Wise et al., 2004), and a time course of neuronal activity (⁠$neu(t),$ varied in magnitude (m), frequency (f), and phase shift (p) relative to $CO2(t)$⁠):

In these simulations (Fig. 7), CVR_BOLD was calculated similarly to as described in Eq. 11, with $CO2(t)$ as the input regressor. For simplicity, a CVR_BOLD of 1 indicates that quantification solely reflects the contribution of the CO₂ time course.

The simulations indicate that CVR_BOLD depends on the magnitude, phase, and frequency properties of resting-state neuronal fluctuations relative to the CO₂-based fluctuations. The influence of resting-state neuronal fluctuations on CVR_BOLD quantification increases with increased resting-state neuronal fluctuation magnitude (m). Calculated CVR_BOLD can either increase (>1) or decrease (<1) depending on the phase/frequency properties of the contaminant neuronal fluctuations. In general, the neuronal fluctuations most dramatically affect CVR_BOLD estimation from 0.015 to 0.075 Hz, a typical range of activity in resting-state networks (often termed the slow-4 and slow-5 bandwidths) and similar in frequency to the simulated CO₂ fluctuations (Buzsáki & Draguhn, 2004; Penttonen & Buzsáki, 2003; Wise et al., 2004). These simulations stress the value of incorporating breath-hold modulations into the resting-state scan, as has been performed in previous work, to introduce larger CO₂ changes and minimize the relative contribution (i.e., m) of neuronal activity (Liu et al., 2020; Stickland et al., 2021, 2022).

Of all clinical cases, steno-occlusive disease is the most widely studied using CVR_BOLD (Sleight et al., 2021). Using simulations, we explored the relationship between CVR_BOLD and arterial stenosis in both the presence and absence of collateral circulation (Sobczyk et al., 2020) to determine whether and to what extent CVR_BOLD is representative of CVR impairment.

In summary (see Section 2.1 of the Supplementary Materials for details), we simulated a tissue voxel supplied by two independent arterial vessels—Vessel A (with variable stenosis; x-axis) and Vessel B (with variable collateral flow; blue arrow) (Fig. 8). The CO₂ time course input to the simulated tissue is thus a summation of the dispersed CO₂ time course from Vessel A and the undispersed CO₂ time course from Vessel B. The tissue was simulated as previously described (see Section 2.2), with one notable modification: to account for autoregulatory exhaustion, tissue vessels were only simulated to respond to the CO₂ stimulus if the cumulative upstream baseline flow (i.e., summated flow from vessels A and B) was greater than 50% of flow in non-stenotic Vessel A (Lassen, 1959).

Many studies report low CVR_BOLD (Duffin et al., 2018; Hartkamp et al., 2017, 2018; Venkatraghavan et al., 2018) and long delay/dispersion (Duffin et al., 2015; Hartkamp et al., 2012; Waddle et al., 2020) in the ipsilateral cortex of steno-occlusive disease. The simulations (Fig. 8B), prior to reaching autoregulatory capacity (i.e., before red arrow), indicate that CVR_BOLD decreases due to stenosis-mediated dispersion of the input CO₂ bolus in the absence of any changes to the tissue vasculature’s physiological CVR. For example, CVR_BOLD reduces by roughly 50% (i.e., from 0.057 s^-1/mmHg to 0.028 s^-1/mmHg) when the simulated upstream flow is reduced by 40%, even though the tissueCVR_Factoris unchanged.