Abstract
Cerebrovascular reactivity (CVR) imaging is used to assess the vasodilatory capacity of cerebral blood vessels. While blood flow (CVRCBF), blood velocity (CVRv), and preferably blood volume changes (CVRCBV) are used to represent physiological CVR, quantifying these measures is fraught with acquisition challenges in humans. Consequently, blood oxygenation level-dependent (BOLD)-MRI CVR (CVRBOLD) is the most widely used MRI-based CVR method, even though it arguably provides the most indirect estimation of CVR. In this paper, we sought to holistically address the quantitative capacity and shortcomings of CVRBOLD. To do so, we developed a CVRBOLD simulation framework and, together with data from the CVRBOLD literature, addressed whether and to what extent CVRBOLD accurately reflects CVR, and with which parameters CVRBOLD varies most. In short, we show the following: CVRBOLD does not necessarily correspond to physiological measures of CVR and depends on physiological (e.g., hematocrit) and acquisition (e.g., field strength) parameters; CVRBOLD is dependent on the stimulus protocol (e.g., breath-holding vs. controlled hypercapnia) chosen to elicit a vasoactive response; resting-state CVRBOLD does not necessarily reflect breath-hold CVRBOLD, likely due to confounding neuronal activity; in stenotic disease and steal physiology, CVRBOLD results from a combination of factors which do not necessarily reflect the underlying CVR. We are confident that this work will provide researchers and clinicians with invaluable insights and advance the field of cerebrovascular imaging by enabling more accurate quantification of CVR in both health and disease.
1 Background
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 CO2, and to an extent, hypoxia (Kety & Schmidt, 1948; Mardimae et al., 2012). Specifically, it is well known that increased partial pressure of CO2 (PCO2) 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 CO2 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 CO2 (i.e., ) in the vessel(s) of interest. This index () 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, 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 and , 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 (CVRBOLD), which is used to measure T2*- or T2-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, CVRBOLD 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 CO2 stimulus protocols, which either require the manipulation of or the harnessing of endogenous fluctuations. While earlier studies focused predominantly on the former approach (e.g., breath-holding, rebreathing, and prospective end-tidal targeting), resting-state signal-based CVRBOLD 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 CVRBOLD simulation framework and interrogated the CVRBOLD literature to address whether and to what extent CVRBOLD accurately reflects ground-truth CVR. First, we investigated the relationship between quantified CVRBOLD 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 CVRBOLD relative to one another. Finally, we simulated example cases of vascular pathology to assess the degree to which CVRBOLD reflects ground-truth CVR in disease.
2 Theory
2.1 Physiological CVR theory and simulations
As previously described, the ground-truth (i.e., most direct) definition of CVR is . 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.
Variation among physiological CVR measures. (A) Simulated CVRs (%/mmHg) as a function of . Here, CVR measures are normalized to an expected of ~9.25 mmHg associated with a of 40% (Grüne et al., 2015; Poulin et al., 1996; Ramsay et al., 1993; Sato et al., 2012; Tancredi & Hoge, 2013); α = 0.29 (Ito et al., 2003). (B) Simulated CVRs as a function of (to represent variations in healthy vessels within the vascular network); = 1.
Variation among physiological CVR measures. (A) Simulated CVRs (%/mmHg) as a function of . Here, CVR measures are normalized to an expected of ~9.25 mmHg associated with a of 40% (Grüne et al., 2015; Poulin et al., 1996; Ramsay et al., 1993; Sato et al., 2012; Tancredi & Hoge, 2013); α = 0.29 (Ito et al., 2003). (B) Simulated CVRs as a function of (to represent variations in healthy vessels within the vascular network); = 1.
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 T1 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 CO2 ( (mmHg)), the corresponding global CBF percent change , and a CVR heterogeneity factor (). The is a novel term used in this work to represent spatial heterogeneity in physiological CVR across the cerebral vasculature. The represents the factor increase (> 1) or decrease (< 1) in percent flow change for a given simulated voxel relative to the global in the brain (i.e., if = 2, simulated voxel corresponds to double the , normalized to ). Thus, is effectively a proxy to voxel-specific , scaled by and .
To quantify (Eq. 2.2) in terms of , Grubb’s relationship is used (Grubb et al., 1974). then depends on , the , and Grubb’s exponent () which governs the relationship between CBF and CBV changes as .
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 , 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 (Eq. 2.3) in terms of , 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 or α. For Figure 1A, the simulated vessels for each 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 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 is ~0.4%/mmHg (when normalized to the expected during breath-holding) (Donahue et al., 2009; Lu et al., 2003). While is similarly expected to be lower than , 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 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 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, , , and will be discordant when comparing vessels with different values (i.e., a voxel with a higher value (arterially dominated) will yield a higher and lower than in a voxel with a lower value (venous dominated), even though may be identical across these vessels, if for example, they are all within the same vascular path). Therefore, and will likely not adequately represent regional variations in the ground-truth (i.e., ) 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, CVRBOLD, an index most widely utilized for CVR measurements, albeit beholden to the confounds of BOLD-MRI, will be the focus of our discussion.
2.2 CVRBOLD theory and simulations
BOLD-MRI is used to measure T2- or T2*-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).
CVRBOLD is typically measured as a percent change in T2*-weighted signal normalized to (Fisher et al., 2018; Liu et al., 2019):
Here, is the baseline BOLD signal (prior to the CO2 challenge) and is the maximum BOLD signal in response to the CO2 challenge. CVRBOLD is assumed to reflect CVR, particularly , given that increased CBF leads to a decrease in deoxyhemoglobin (dOHb) in the capillaries and veins, and thus, an increase in T2*-weighted MRI signal (as summarized in Liu et al., 2019 and Uludağ et al., 2009). However, to link with CVRBOLD, the contribution of various physiological and physical variables must be considered and assessed.
2.2.1 Simulations: Signal model overview
Using a simulation framework, we sought to determine the influence that various parameters (e.g., echo time, field strength, blood volume) have on CVRBOLD 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).
T2*-weighted MRI signal () is a summation of extravascular (EES) and intravascular (IVS) signal contributions, and , 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 ( (%)), each vessel’s volume fraction (, reflecting the fraction of blood volume attributable to a given vessel type), and each vessel’s hypercapnia-induced decimal change in blood volume () (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, in our simulations does not consider spin density effects, but this can easily be incorporated in the future by modifying ϕ. In addition, is normalized such that at an echo time of 0 ms, is equal to 1.
The signal components are separable as follows:
and (s-1) are the extravascular and intravascular magnetic relaxation rates induced by a susceptibility agent, respectively. and (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).
2.2.2 Simulations: Baseline tissue/vascular properties
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 CVRBOLD quantification.
The extravascular and intravascular relaxation rates across field strength found in the literature.
. | [s-1] . | [s-1] . | [s-1] . | [s-1] . |
---|---|---|---|---|
1.5 T (Hct=44%) | 15.38 | 7.5 | Capillary: | |
Arteriole/Venule: | ||||
3 T (Hct=44%) | 18.03 | 20.7 | Capillary: | |
Arteriole/Venule: | ||||
4.7 T (Hct=44%) | 27.35 | 42.5 | Capillary: | |
Arteriole/Venule: | ||||
7 T (Hct=44%) | 32.6 | 116 | Capillary: | |
Arteriole/Venule: | ||||
3 T (Hct=21%) | 18.03 | 18 | Capillary: | |
Arteriole/Venule: |
. | [s-1] . | [s-1] . | [s-1] . | [s-1] . |
---|---|---|---|---|
1.5 T (Hct=44%) | 15.38 | 7.5 | Capillary: | |
Arteriole/Venule: | ||||
3 T (Hct=44%) | 18.03 | 20.7 | Capillary: | |
Arteriole/Venule: | ||||
4.7 T (Hct=44%) | 27.35 | 42.5 | Capillary: | |
Arteriole/Venule: | ||||
7 T (Hct=44%) | 32.6 | 116 | Capillary: | |
Arteriole/Venule: | ||||
3 T (Hct=21%) | 18.03 | 18 | Capillary: | |
Arteriole/Venule: |
Relaxation rate values and equations are from Zhao et al. (2007) and Uludağ et al. (2009) at multiple field strengths. Hct (hematocrit) represents the percentage of blood containing red blood cells (see Section 4.4 for Hct simulations). is the decimal baseline oxygen saturation of hemoglobin (see Table 2 for values in each vascular compartment). is the change in oxygen saturation of hemoglobin (see Section 2.2.3 for modeling). TEs used for 1.5 T, 3 T, 4.7 T, and 7 T were 40 ms, 30 ms, 25 ms, and 20 ms, respectively.
Summary of simulated tissue/vascular properties.
. | Tissue . | ||
---|---|---|---|
. | 4 . | ||
(%) . | Arteriole . | Capillary . | Venule . |
0.2 | 0.4 | 0.4 | |
0.95 | 0.81 | 0.67 |
. | Tissue . | ||
---|---|---|---|
. | 4 . | ||
(%) . | Arteriole . | Capillary . | Venule . |
0.2 | 0.4 | 0.4 | |
0.95 | 0.81 | 0.67 |
The values in this table can be used for any arbitrary set of physiological assumptions (e.g., here tissue represents gray matter (GM), but can be changed to 2% to simulate white matter (WM) as well). Note that changing baseline affects scaling but not the qualitative results.
2.2.3 Simulations: Modeling oxygenation change
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 () can be estimated from the well-known cerebral metabolic rate of oxygen (CMRO2) and CBF relationship (Mintun et al., 1984)—note that CMRO2 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. depends on the simulated vessel. for vein/venule () depends on , the oxygenation extraction (baseline arteriole oxygenation ()—baseline venule oxygenation ()), and the in simulated tissue. To calculate , we first set to 9.25 mmHg based on typical maximal changes during -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 in response to a of 9.25 mmHg ( of ~40%). As previously discussed, is a novel term used in this work to represent spatial heterogeneity in across the cerebral vasculature (i.e., simulated voxels).
for artery/arteriole () 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).
for capillary () 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 () for a given vessel can then be described as the summation of its and the oxygenation change resulting from hypoxia ().
is assumed to be 0 for all simulations, except for breath-holding (see Section 5).
2.2.4 Simulations: Modeling vasodilation
As previously mentioned (Eq. 4), hypercapnia-induced vasodilation () must be included in the CVR modeling (Kety & Schmidt, 1948; Liu et al., 2019). depends on the , , 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 is relative (e.g., if = 1.08 (i.e., an 8% change) and = 4% (typical values for gray matter), CBV would increase to 4.32%).
2.2.5 Simulations: Calculation of
According to a recent systematic review (Sleight et al., 2021), the vast majority of papers report CVRBOLD as percent signal change normalized to (), described in Eq. 3, whereas only 0.4% of studies report CVRBOLD as a relaxation rate change normalized to ():
As a reminder, is the baseline signal (when is set to 0), is the maximum signal (e.g., when is set to 40%), and is the corresponding change in relaxation rate. The preponderance of relative to is very surprising given the well-known relationship between echo time and T2* signal change (Eqs. 5.1 and 5.2). That is, higher TE yields higher calculated values of (e.g., see figure 1 in Havlicek et al., 2017), whereas should be (almost) TE independent. As previously mentioned, is set to 9.25 mmHg for these simulations.
3 CVRBOLD Dependency on Acquisition Parameters
3.1 Echo time
In addition to investigating the effect of TE on CVRBOLD using our simulations, we surveyed gray matter (GM) and whole brain (WB) values reported in the literature across a range of TEs at 3 T (Fig. 2).
Relationship between echo time and estimation. (A) Literature: Table and plot of reported values (%/mmHg) in average GM and WB across a range of echo times (ms). (B) Simulations: Relationship between and echo time. (C) Simulations: Relationship between and echo time. Note that the y-axis range is much smaller relative to B. CBV0 was varied from 2% to 32%. Note that scaled CVRBOLD represents the percent change of CVRBOLD in relation to CVRBOLD at 20 ms (i.e., ).
Relationship between echo time and estimation. (A) Literature: Table and plot of reported values (%/mmHg) in average GM and WB across a range of echo times (ms). (B) Simulations: Relationship between and echo time. (C) Simulations: Relationship between and echo time. Note that the y-axis range is much smaller relative to B. CBV0 was varied from 2% to 32%. Note that scaled CVRBOLD represents the percent change of CVRBOLD in relation to CVRBOLD at 20 ms (i.e., ).
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 (Fig. 2A), in agreement with our simulations (Fig. 2B). In fact, Triantafyllou et al. specifically investigated the dependency of 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 across a range of TEs in the literature, our simulations show that the TE-dependency is largely eliminated when CVR is calculated as (i.e., values differ by only ~2–4% between TE = 20 ms and TE = 40 ms) (Fig. 2C).
3.2 Field strength
An increased contrast-to-noise ratio makes a compelling case for the implementation of increased magnetic field strength in CVRBOLD studies; however, the relationship between field strength and both and must be considered.
Triantafyllou et al. investigated 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 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 CVRBOLD, ASL-derived CVR (CVRASL), 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 and (Fig. 3).
Relationship between field strength and CVRBOLD estimation. (A) Literature: Table and plot of reported values (%/mmHg) in average GM across a range of field strengths. (B) Simulations: Relationship between and field strength. (C) Simulations: Relationship between and field strength. CBV0 was varied from 2 to 32%. Note that scaled CVRBOLD represents the percent change of CVRBOLD in relation to CVRBOLD at 1.5 T (i.e., ).
Relationship between field strength and CVRBOLD estimation. (A) Literature: Table and plot of reported values (%/mmHg) in average GM across a range of field strengths. (B) Simulations: Relationship between and field strength. (C) Simulations: Relationship between and field strength. CBV0 was varied from 2 to 32%. Note that scaled CVRBOLD represents the percent change of CVRBOLD in relation to CVRBOLD at 1.5 T (i.e., ).
GM 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 (R2 = 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 CVRBOLD 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 and as a function of field strength—in theory, this would likely result in regional scaling differences between CVRBOLD 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 CVRBOLD measures. Upon fitting the and simulations, we found that logarithmic-radical fits best corrected for (R2 = 0.999) and (R2 = 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 CVRBOLD values to ensure better inter-field strength comparability. Ideally, this correction factor would be derived and validated empirically by collecting CVRBOLD 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 quantifyCVRBOLDas for the remainder of this work. Although not explicitly shown, results are qualitatively identical to in the forthcoming relationships.
4 CVRBOLD Dependency on Physiological Parameters
The quantitative relationships between CVRBOLD 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 CVRBOLD on a comprehensive set of physiological parameters.
4.1 Baseline blood volume
One parameter that is expected to influence CVRBOLD is the underlying baseline CBV (CBV0), in agreement with the literature (Blockley et al., 2013; Buxton et al., 2004; Davis et al., 1998; Yablonskiy & Haacke 1994), where the relationship between CBV0 and the percent BOLD change (and/or ) 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 CBV0, 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 CBV0 results in a doubling of the measured BOLD-MRI signal change, and thus, CVRBOLD. Consequently, comparing CVRBOLD,GM with CVRBOLD,WM, for example, is heavily biased by the roughly doubled CBV0 in GM (CBV0 ~ 4%) relative to WM (CBV0 ~ 2%). Therefore, making the claim that “CVRBOLD reflects ” is just as valid as claiming that “CVRBOLD reflects CBV0,” which poses a considerable limitation on CVRBOLD in its current state. To clarify, if CVRBOLD reflected CVRCBV (that is, the change in CBV), this would be a benefit, not a limitation. However, CVRBOLD scales linearly with baselineCBV(CBV0). 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), CVRBOLD,GM has empirically been found to be roughly double CVRBOLD,WM in the literature, while CVRASL,GM has been found to be roughly the same as CVRASL,WM (Taneja et al., 2020; Xu et al., 2024). Since DSC-MRI calculated CBV0 (Ostergaard, 2005; Schulman et al., 2023) and CVRBOLD both scale roughly linearly with CBV0, dividing CVRBOLD by CBV0 (for example, measured separately using DSC) is a promising solution for removing the confound of CBV0 from measures of CVRBOLD. 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 CVRASL which is assumed to be independent of CBV0 (Petersen et al., 2006; Taneja et al., 2020).
4.2 Arteriolar blood volume
The observed tissue relaxation rate change in CVRBOLD 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 CBV0 (CBVA (%)) on the measured CVRBOLD (Fig. 4A).
Simulated relationship between CVRBOLD and blood/tissue parameters. Scaled CVRBOLD represents the percent change of CVRBOLD in relation to CVRBOLD at the first datum point for subplots A–B (e.g., in (A): ). (A) Relationship between relative CBVA (percentage of CBV0 that is arteriolar (like arteriolar Vi in Section 2.2.2, but as a percentage)) and estimated CVRBOLD. The was varied from 0.5 to 2.5 for subplot A. (B) Relationship between the baseline arteriolar blood oxygenation () and estimated CVRBOLD. CBV0 was varied from 2 to 32%. (C) Relationship between the Hct (%) and estimated CVRBOLD. Black represents control (Hct = 44%) and cyan represents anemia (Hct = 21%). Scaled CVRBOLD represents the percentage difference of CVRBOLD in relation to CVRBOLD for control, when CVRFactor = 1.
Simulated relationship between CVRBOLD and blood/tissue parameters. Scaled CVRBOLD represents the percent change of CVRBOLD in relation to CVRBOLD at the first datum point for subplots A–B (e.g., in (A): ). (A) Relationship between relative CBVA (percentage of CBV0 that is arteriolar (like arteriolar Vi in Section 2.2.2, but as a percentage)) and estimated CVRBOLD. The was varied from 0.5 to 2.5 for subplot A. (B) Relationship between the baseline arteriolar blood oxygenation () and estimated CVRBOLD. CBV0 was varied from 2 to 32%. (C) Relationship between the Hct (%) and estimated CVRBOLD. Black represents control (Hct = 44%) and cyan represents anemia (Hct = 21%). Scaled CVRBOLD represents the percentage difference of CVRBOLD in relation to CVRBOLD for control, when CVRFactor = 1.
While and are held constant, CVRBOLD 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 CVRBOLD 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, CVRBOLD has a low sensitivity to CVR in arteries/arterioles (see Schulman et al., 2024, which shows how CVRBOLD 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 CVRBOLD.
4.3 Baseline blood oxygenation
While baseline arterial/arteriolar blood oxygenation ( (%)) 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 on CVRBOLD (Fig. 4B). Of note, it is known that reducing beyond roughly 80% (i.e., ~50 mmHg) results in a hypoxia-driven cerebrovascular response (Mardimae et al. 2012); therefore, we constrained the parameter to avoid simulating this additional hypoxia-mediated increase in CBF.
The relationship between and CVRBOLD is heavily dependent on the CBV0 of the simulated voxel (i.e., CVRBOLD magnitude increases non-linearly with in tissue (CBV0 ~2–8%) but decreases non-linearly with in highly vascular voxels (e.g., CBV0 of 32%)). The observation for the highly vascular voxels is in agreement with our previous DSC MRI work, which found that induced by susceptibility agent was negatively correlated with 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 CVRBOLD increase as a function of in tissue is owed to a larger overall extravascular relaxation rate change magnitude at higher 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 (refer to Schulman et al., 2024). In general, these simulations reveal a complex relationship between the measured CVRBOLD and which is important for consideration in hypoxemic patients. In future work, incorporating the measured (estimated from a patient’s SpO2) into tissue CVRBOLD modeling may bypass this dependency in patients with abnormally low values.
4.4 Hematocrit
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 CVRBOLD (Kosinski et al., 2017), it is unclear to what extent Hct attenuates CVRBOLD when and are unchanged. To examine this relationship, we simulated (see Section 2.2.1 and Table 1) CVRBOLD 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., = 1), the quantified CVRBOLD is reduced in simulated anemia by ~40%, relative to normal tissue—this appears regardless of whether the is low (0.5) or high (2). These results have important implications when assessing CVRBOLD in anemic patients, and even when comparing results across sex—namely, while reduced might partially account for the observed CVRBOLD 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 CVRBOLD increases as a function of Hct even when physiological CVR (i.e., or ) is unchanged. Upon validating this dependency experimentally (i.e., observing CVRBOLD 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 CVRBOLD dependency.
Note that variations in baseline CMRO2 would theoretically confound the relationship between Hct and CVRBOLD as it would alter (i.e., increased baseline CMRO2 would yield a decrease in ). CVRBOLD is dependent on , as we have shown above (Fig. 4B).
5 CVRBOLD Dependency on Stimulus Design
5.1 Controlled hypercapnia versus breath-holding
Breath-holding is a well-accepted method in CVRBOLD 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 CVRBOLD values.
Using hypoxia and hypercapnia values reported in previous works (refer to Supplementary Fig. 1 and Supplementary Table 1), we simulated CVRBOLD 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 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 CVRBOLD values to (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.
Simulated relationship between CVRBOLD estimation and breath-hold duration. CVRBOLD normalized (red; s-1/mmHg) or not normalized (blue; s-1) to (Supplementary Table 1 and Supplementary Fig. 1). For quantification, see Eq. 9.2 and Supplementary Table 1 for corresponding values. (A) Simulated without hypoxia (to mimic controlled hypercapnia). (B) Simulated with hypoxia (to mimic breath-holding).
Simulated relationship between CVRBOLD estimation and breath-hold duration. CVRBOLD normalized (red; s-1/mmHg) or not normalized (blue; s-1) to (Supplementary Table 1 and Supplementary Fig. 1). For quantification, see Eq. 9.2 and Supplementary Table 1 for corresponding values. (A) Simulated without hypoxia (to mimic controlled hypercapnia). (B) Simulated with hypoxia (to mimic breath-holding).
In the controlled hypercapnia simulations, -normalized CVRBOLD progressively decreases, though minimally, and -unnormalized CVRBOLD progressively increases with the stimulus magnitude until saturates after changing by ~10.5 mmHg (Sasse et al., 1996) (Fig. 5A). The mild decrease in normalized CVRBOLD is initially perplexing, especially given that mildly increases as a function of (Supplementary Table 2). This becomes clearer when recognizing that BOLD begins to saturate with progressively larger increases in and , 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, -normalized CVRBOLD decreases substantially, and -unnormalized CVRBOLD 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 -unnormalized CVRBOLD 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 CVRBOLD measurement (Bright & Murphy, 2013); here, CVRBOLD 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 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 values measured from Bright and Murphy (2013) into our model (instead of those from Sasse et al., 1996), we obtain -normalized CVRBOLD 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 CVRBOLD (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 CVRBOLD 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 CVRBOLD experiments. However, the simulation results suggest that those seeking to limit the hypoxic effect (i.e., those seeking breath-hold CVRBOLD values in agreement with controlled-hypercapnia CVRBOLD values) should aim for a breath-hold duration of ~10–15 s—this range is anyway expected to be more translatable to patients.
5.2 Breath-holding versus resting-state
Over the past two decades, researchers have claimed that CVRBOLD 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 fluctuations (Wise et al., 2004). In general, resting-state CVRBOLD 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 CVRBOLD 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 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).
5.2.1 Experimental data
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 CVRBOLD (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 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 (). 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 was fit to each tissue voxel by first shifting the forward in 0.5 s steps, up to a maximum of 7 s, and the correlation between the 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 was then fit to the time course in each voxel (), using the voxel’s delay as an input parameter, assuming the following model:
Here, is the scaling factor that minimizes least squared error between the tissue and relaxation rate time courses. In terms of the resulting maps, is effectively a proxy to a regression-based estimation of (the absolute value of will be different from as values are not normalized to , but comparing absolute values is not the focus of our analysis). Note that is not a proxy to steady-state, which requires the additional modeling of a hemodynamic response function (HRF) (Poublanc et al., 2015). Finally, was normalized to the average GM to obtain the relative CVR maps (). represents the vertical shift needed to account for baseline MRI signal differences and represents the residual error.
In Figure 6, the subject-averaged maps (which are normalized to average GM) show agreement between field strengths in most regions (i.e., patterns of 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 maps (Fig. 6B). In comparison with breath-holding, resting-state yields higher in the visual cortex, pre-central gyrus, and post-central gyrus, but lower in the subcortical GM (e.g., putamen).
Differences between resting-state and breath-hold . (A) and input regressor correlation maps in MNI152 2 mm anatomical space for subject-averaged breath-hold and resting-state data at 3 T (n = 9) and 7 T (n = 9). (B) Statistical comparison (paired t-test, α = 0.05) between breath-hold and resting-state data at 7 T (see Supplementary Fig. 3 for 3 T results). Only significant t-scores are displayed; RS > BH corresponds with t > 0.
Differences between resting-state and breath-hold . (A) and input regressor correlation maps in MNI152 2 mm anatomical space for subject-averaged breath-hold and resting-state data at 3 T (n = 9) and 7 T (n = 9). (B) Statistical comparison (paired t-test, α = 0.05) between breath-hold and resting-state data at 7 T (see Supplementary Fig. 3 for 3 T results). Only significant t-scores are displayed; RS > BH corresponds with t > 0.
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 CO2-induced relaxation fluctuations (Wise et al., 2004), are the most likely explanation for these differences between breath-hold and resting-state 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 . 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 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 CVRBOLD. Although voxel-wise resting-state CVRBOLD 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 CO2-based CVRBOLD maps (i.e., regional correlations below 0.7), with a similar observation of increased resting-state CVRBOLD in the visual cortex relative to traditional CVRBOLD (Liu et al., 2021). Although not calculated here, we also expect differences between breath-hold and resting-state absolute values due to factors such as arterial compliance (Prokopiou et al., 2019).
As a final note, breath-holding-based CVRBOLD 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 CVRBOLD is similarly confounded by neuronal activity in certain brain regions.
5.2.2 Simulations
To better understand the influence of concomitant neuronal activity on CVRBOLD 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 CO2 time course (, simple sinusoid with a frequency of 0.033 Hz Wise et al., 2004), and a time course of neuronal activity ( varied in magnitude (m), frequency (f), and phase shift (p) relative to ):
In these simulations (Fig. 7), CVRBOLD was calculated similarly to as described in Eq. 11, with as the input regressor. For simplicity, a CVRBOLD of 1 indicates that quantification solely reflects the contribution of the CO2 time course.
Simulated effect of resting-state neuronal activity on CVRBOLD quantification. An input CO2-based time course and tissue time course (summation of the input time course and a time course of neuronal activity (varied in frequency (y-axis), phase (x-axis), and magnitude (m) relative to the input time course)) were simulated. CVRBOLD was calculated as described in the experimental data (i.e., input time course was parameterized with delay and CVRBOLD (effectively a magnitude scaling factor), and least squares-based optimization was then conducted between this parameterized time course and the tissue time course).
Simulated effect of resting-state neuronal activity on CVRBOLD quantification. An input CO2-based time course and tissue time course (summation of the input time course and a time course of neuronal activity (varied in frequency (y-axis), phase (x-axis), and magnitude (m) relative to the input time course)) were simulated. CVRBOLD was calculated as described in the experimental data (i.e., input time course was parameterized with delay and CVRBOLD (effectively a magnitude scaling factor), and least squares-based optimization was then conducted between this parameterized time course and the tissue time course).
The simulations indicate that CVRBOLD depends on the magnitude, phase, and frequency properties of resting-state neuronal fluctuations relative to the CO2-based fluctuations. The influence of resting-state neuronal fluctuations on CVRBOLD quantification increases with increased resting-state neuronal fluctuation magnitude (m). Calculated CVRBOLD 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 CVRBOLD 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 CO2 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 CO2 changes and minimize the relative contribution (i.e., m) of neuronal activity (Liu et al., 2020; Stickland et al., 2021, 2022).
6 Disease Modeling and CVRBOLD
6.1 Vascular occlusion and collateral flow
Of all clinical cases, steno-occlusive disease is the most widely studied using CVRBOLD (Sleight et al., 2021). Using simulations, we explored the relationship between CVRBOLD and arterial stenosis in both the presence and absence of collateral circulation (Sobczyk et al., 2020) to determine whether and to what extent CVRBOLD 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 CO2 time course input to the simulated tissue is thus a summation of the dispersed CO2 time course from Vessel A and the undispersed CO2 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 CO2 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).
Simulated relationship between CVRBOLD estimation, vascular occlusion, and collateral flow. (A) Schematic of stenotic (Vessel A) and collateral (Vessel B) vessels supplying the imaged tissue voxel. (B) CVRBOLD (s-1/mmHg) for simulated tissue (downstream of Vessels A and B) that was capable of autoregulation 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. Blue arrow indicates increasing collateral flow contribution (i.e., Vessel B flow simulated from 0 to 100% of the flow supplied by a non-stenotic Vessel A). Red arrow indicates when autoregulatory capacity is reached.
Simulated relationship between CVRBOLD estimation, vascular occlusion, and collateral flow. (A) Schematic of stenotic (Vessel A) and collateral (Vessel B) vessels supplying the imaged tissue voxel. (B) CVRBOLD (s-1/mmHg) for simulated tissue (downstream of Vessels A and B) that was capable of autoregulation 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. Blue arrow indicates increasing collateral flow contribution (i.e., Vessel B flow simulated from 0 to 100% of the flow supplied by a non-stenotic Vessel A). Red arrow indicates when autoregulatory capacity is reached.
Many studies report low CVRBOLD (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 CVRBOLD decreases due to stenosis-mediated dispersion of the input CO2 bolus in the absence of any changes to the tissue vasculature’s physiological CVR. For example, CVRBOLD 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 tissueCVRFactoris unchanged.
Dispersion effects are not only relevant in disease, but also when comparing healthy tissue with different dispersion/transit properties, such as the GM and WM. As demonstrated in DSC-MRI experiments, WM is associated with more dispersion of the input bolus (i.e., from artery to tissue) than GM (Ibaraki et al., 2007). From this, it can be expected that the observed CVRBOLD will be reduced in WM using a traditional CVR analysis, even though this dispersion-mediated reduction in bolus peak does not reflect a reduction in the tissue’s ability to vasodilate (i.e., the tissue is simply seeing a reduced CO2 peak due to dispersion and will yield a smaller response accordingly). Given this CVRBOLD dependency, it is recommended that researchers calculate the steady-state CVRBOLD by implementing a model-based hemodynamic response function (Poublanc et al., 2015) which minimizes the dispersion confound, or calculate CVRBOLD as the integral of the CVRBOLD time course—a potential model-free, computationally inexpensive alternative, similar to the CBV0 calculation in DSC-MRI (Bjørnerud & Emblem, 2010; Meier & Zierler, 1954; Ostergaard et al., 1996).
Once the tissue arterioles reach autoregulatory capacity (Fig. 8B, red arrow), CVRBOLD is 0. In fact, although not shown in this simulation, the co-occurrence of steal physiology could even yield a negative CVRBOLD measurement when the tissue arterioles have reached autoregulatory capacity (Duffin et al., 2017; Fisher & Mikulis, 2021; Sobczyk et al., 2014).
The degree of stenosis-mediated CVRBOLD reduction is heavily modulated by the degree of collateral flow (i.e., in a scenario where collateral flow replaces the entirety of flow from Vessel A, CVRBOLD would be in agreement with CVRBOLD in the absence of both occlusion and collateral flow). This agrees with the literature, where collateral flow has been posited as an explanation for the counterintuitive positive CVRBOLD observed in vascular territories affected by steno-occlusion (Sobczyk et al., 2020). Of note, increased collateral flow helps prevent the simulated tissue from reaching autoregulatory capacity (see Fig. 8B, where vasodilatory exhaustion does not occur with collateral flow greater than 50% of flow relative to non-stenotic Vessel A).
In sum, the simulations indicate that CVRBOLD may reflect physiological CVR impairment when tissue vessels have reached autoregulatory capacity, but not in the case where there is a remaining vasodilatory reserve due to stenosis-mediated bolus dispersion; this limitation may be partially attenuated by using more advanced CVRBOLD analysis strategies that model the effects of dispersion (Bjørnerud & Emblem, 2010; Meier & Zierler, 1954; Ostergaard et al., 1996; Poublanc et al., 2015).
6.2 Steal physiology
Steal physiology describes the phenomenon in which vasodilatory capacity is exhausted in one vessel, such that parallel, unaffected vessels “steal” blood flow away from the exhausted vessel when presented with a vasodilatory stimulus (Brawley, 1968; Symon, 1968). CVRBOLD has often been used to detect steal physiology in the brain, where negative CVRBOLD values are often assumed to indicate regions where the blood vessels see a reduction in flow in the presence of a vasodilatory stimulus (Fisher & Mikulis, 2021; Sobczyk et al., 2014). It is not entirely clear to what extent CVRBOLD-measured steal physiology reflects ground-truth CVR impairments. To address this, we simulated steal physiology as a circuit with parallel flow and resistance (Fig. 9A and Supplementary Fig. 7). Note that these simulations are completely independent from those described in Section 2.2.
Steal physiology: Simulated relationship between vascular reactivity and flow change. (A) Schematic of steal physiology. Upstream vasculature unable to dilate due to stenosis. Vessel 3 vasodilatory capacity depends on the level of downstream stenosis. (B) Percent flow change in Vessel 3 (in the presence of a vasodilative agent) as a function of vascular reactivity in Vessel 3 (y-axis) and vascular reactivity in Vessel 2 (x-axis). Four different types of CVR responses are indicated here (see text).
Steal physiology: Simulated relationship between vascular reactivity and flow change. (A) Schematic of steal physiology. Upstream vasculature unable to dilate due to stenosis. Vessel 3 vasodilatory capacity depends on the level of downstream stenosis. (B) Percent flow change in Vessel 3 (in the presence of a vasodilative agent) as a function of vascular reactivity in Vessel 3 (y-axis) and vascular reactivity in Vessel 2 (x-axis). Four different types of CVR responses are indicated here (see text).
In summary (see Section 2.2 of the Supplementary Materials for details), this simulation was set up as a parallel flow circuit (Fig. 9A and Supplementary Fig. 7). Baseline resistances in Vessels 2 and 3 were set to 1; baseline resistances in Vessels 1 and 4, which represents a summation of vasculature prior to and after the parallel circuit, respectively, were set to 1/3. As “voltage” was not found to affect the quantitative values in this simulation, it was set to 1. Using Kirchoff’s laws, the resulting flow through Vessel 3 (and thus, the supplied tissue vasculature) was calculated iteratively by reducing resistance (from 100 to 0%, in 1% steps) in Vessels 2 and 3 (x- and y-axis, respectively, in Fig, 9B). We assume that Vessel 4 (i.e., venous vasculature) sees a negligible resistance change relative to the other vessels (Grubb et al., 1974; Ito et al., 2003) and Vessel 1, due to simulated occlusion, also sees a negligible change—for the scenario where Vessel 1 is capable of vasodilation, see Supplementary Figure 8.
In Figure 9B, note that positive percent flow change corresponds to positive in tissue supplied by Vessel 3, whereas negative percent flow change corresponds to negative in tissue supplied by Vessel 3. For simplicity here, the “sign” of can be assumed to be the same sign for CVRBOLD; therefore, flow change can be thought of as a proxy to CVRBOLD.
Reducing resistance in Vessel 2 (via CO2-mediated dilation) while maintaining constant resistance in Vessel 3 (i.e., due to exhausted vasodilatory reserve) yields a negative flow change (i.e., steal physiology) in Vessel 3, and thus, negative CVRBOLD; however, this is not necessarily the case in a scenario where upstream vasculature is capable of vasodilation (Supplementary Fig. 8A). In particular, flow may be relatively unchanged in Vessel 3, even up to a ~40% resistance drop in Vessel 2 (Supplementary Fig. 8B), which is a typical resistance change observed during hypercapnia (Aslan et al., 2010; Blockley et al., 2013; Bulte et al., 2012; De Vis et al., 2015; Gauthier et al., 2013; Sasse et al., 1996; Vestergaard & Larsson, 2019). Thus, any observed negative CVRBOLD in this scenario would not be a result of steal physiology but may instead result from a low contrast-to-noise ratio, vasodilation-mediated negative signal change (Thomas et al., 2013), or errors in analysis (such as not performing lag optimization; Stickland et al., 2021). However, when the upstream vasculature is entirely incapable of vasodilation (Fig. 9B), steal physiology is more likely to be observed, as any resistance decrease in Vessel 2 will yield a negative flow change in Vessel 3, assuming Vessel 3 is incapable of vasodilation.
Interestingly, a steal-mediated decrease in flow may accompany an increase in blood volume for a vessel with some remaining vasodilatory reserve. For example, a 10% resistance decrease in Vessel 3 (i.e., percent blood volume increase of ~2.5%) yields a negative flow change (i.e., negative CVRBOLD) when Vessel 2 resistance decreases by 25% or more (Fig. 9B). Therefore, when studying steal physiology, it is important to identify which representation of physiological CVR (i.e., Eqs. 1.1–1.3) is of greater interest, as agreement between the measures is not guaranteed (e.g., CVRVASO could yield a positive result when CVRASL yields a negative result).
The simulations provide an explanation for all four types of CVR responses which have been documented in previous studies (Duffin et al., 2017; Sobczyk et al., 2014): type A response (i.e., positive flow change with increasing CO2), type B response (i.e., positive flow change with increasing CO2, followed by negative flow change with further CO2 increases), type C response (i.e., negative flow change with increasing CO2), and type D response (i.e., negative flow change with increasing CO2, followed by positive flow change with further CO2 increases) (Fig. 9B).
The simulations demonstrate that CVRBOLD in steal physiology depends on the vasodilatory capacity of multiple vessels, influenced by their respective degrees of stenosis. Importantly, the simulations also demonstrate that CVRBOLD may not always accurately reflect the ground-truth CVR properties of the imaged vessels.
7 Limitations
In this paper, we used a simulation framework (modified from prior fMRI/DSC studies; Kjolby et al., 2006, 2009; Schulman et al., 2022, 2023; Uludağ et al., 2009) to understand how CVRBOLD is influenced by various physiological and physical parameters. Although the model is quite comprehensive, there are a few notable limitations that we would like to address. First, in our simulations, we assume that CMRO2 remains constant during mild/moderate hypercapnic challenges, in agreement with findings from previous work (Blockley et al., 2013; Chen & Pike, 2010; Jain et al., 2011; Vestergaard & Larsson, 2019; Zappe et al., 2008). These studies indicate that CMRO2 does not significantly change during hypercapnia. Fortunately, our modeling framework of CVRBOLD can easily accommodate CMRO2 changes if needed for certain subjects or experimental conditions (in particular, Eq. 6.1). Second, we simulated tissue composed of arteriole, capillary, and venule. However, cerebral tissue may also contain contributions from larger arteries and veins. Future studies focused on understanding the influence of larger vessels on CVRBOLD can still implement our model by incorporating Eqs. 10–11 from Uludağ et al., 2009 for vessel diameters ≥ 200 μm. Third, we assume a ~40% in response to a of 9.25 mmHg, based on previous relationships determined experimentally (Grüne et al., 2015; Poulin et al., 1996; Ramsay et al., 1993; Sato et al., 2012; Tancredi & Hoge, 2013). From study-to-study, the input stimulus may differ substantially from 9.25 mmHg (e.g., in resting-state vs. breath-hold vs. controlled hypercapnia). Therefore, future studies looking to simulate a smaller/larger with our model will need to modify based on the relationships described in the literature (Grüne et al., 2015; Poulin et al., 1996; Ramsay et al., 1993; Sato et al., 2012; Tancredi & Hoge, 2013). Implementing other values will change quantitative values but not alter the qualitative results/findings. Finally, it should be noted that although our simulation framework works as a forward model, using it as an inverse model (i.e., as a replacement for Eq. 3 for the calculation of CVRBOLD) would likely not be possible unless one were to estimate multiple physiological measures (i.e., CBV0, Hct, and Y0) independently from the CVRBOLD acquisition.
In addition to our simulated work, we collected empirical breath-hold and resting-state data to understand whether CVRBOLD is concordant between these methods. There are two notable limitations with respect to the empirical portion of this work. First, data were not available for collection in this study—instead, we used an AIF time course as our input regressor. Scaling signal/relaxation change by is important for obtaining absolute CVRBOLD values. However, it should be noted that in this work, we were not interested in comparing absolute values; rather, we were interested in performing a qualitative voxel-wise comparison between resting-state and breath-hold maps (we do not expect scaling to confound our relative maps). Second, we did not instruct participants to keep their eyes closed or open during the resting-state scan. As has been observed in previous work (Han et al., 2023), this parameter is known to influence resting-state neuronal activity measured with fMRI, which may have contributed to some of the inter-subject variability in our study.
8 Summary
In this work, we sought to understand the differences between several formulations of physiological CVR (i.e., CVRCBV, CVRCBF, and CVRv) and holistically address the accuracy of CVRBOLD. To do so, we developed a CVRBOLD simulation framework and interrogated the CVRBOLD literature to address whether and to what extent CVRBOLD accurately reflects physiological CVR, and with which parameters CVRBOLD varies the most. In summary, we have shown the following:
The various formulations of physiological CVR (CVRCBV, CVRCBF, and CVRv) will not yield the same absolute values, and they do not yield the same relative maps in the presence of regional variations in the composition of blood vessels.
CVRBOLD is linearly dependent on the echo time and non-linearly dependent on the magnetic field strength; quantifying CVRBOLD as CVRBOLD,R largely circumvents the dependency on echo time, while implementing an empirically validated correction factor should minimize the dependency on field strength.
CVRBOLD is dependent on several physiological parameters, including the baseline blood volume, baseline oxygen saturation, hematocrit, and relative contribution of arterial blood; various correction factor-based approaches, which may require acquiring additional physiological data, can be utilized to minimize the influence of these parameters.
Breath-hold-based CVRBOLD may be limited by a competing hypoxic effect; the magnitude of hypoxia increases with breath-hold duration, suggesting that short–moderate breath-hold durations (~15–20 s) will yield a decent balance between increasing the contrast-to-noise ratio and limiting the hypoxic effect.
Resting-state-based CVRBOLD is confounded by the contribution of neuronal activity, possibly creating a discrepancy with CVRBOLD maps obtained using other approaches, such as breath-holding; this contribution can likely be reduced by introducing even basic breathing maneuvers during the scan.
In stenotic disease, CVRBOLD is limited by stenosis-mediated dispersion, which is attenuated by the degree of collateral flow; using a CVRBOLD model which incorporates a localized hemodynamic response function can minimize this limitation.
Negative CVRBOLD cannot always be attributed to steal physiology; likewise, positive CVRBOLD should not always be interpreted as an absence of steal physiology.
In shedding light on these quantitative nuances and providing novel solutions to circumvent some of these problems, we are optimistic that this work will help researchers and physicians more appropriately quantify and interpret CVR in both health and disease.
Data and Code Availability
The dataset used for the current study is available from the corresponding authors upon reasonable request.
Simulation code can be found at the following GitHub link:
Author Contributions
J.B.S.: Conceptualization, Data Curation, Formal Analysis, Investigation, Methodology, Software, Visualization, and Writing. K.U.: Conceptualization, Formal Analysis, Funding Acquisition, Investigation, Methodology, Project Administration, Resources, Supervision, and Writing.
Ethics Statement
The experimental section of this work was approved by the Research Ethics Board of Sungkyunkwan University and all procedures followed the principles expressed in the Declaration of Helsinki. Informed consent was obtained in all 10 healthy volunteers.
Funding
The research conducted in this paper was supported by funding from the Canadian Institutes of Health Research (CIHR) to K.U. and J.S. The study was also supported by the Institute for Basic Science, Suwon, Republic of Korea (IBS-R015-D1) to K.U.
Declaration of Competing Interest
The listed authors have no conflicts of interest to declare.
Acknowledgments
We would like to thank Seong-Gi Kim, Boohee Choi, and Suji Jeong at the Institute for Basic Science for helping us acquire the experimental data.
Supplementary Materials
Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00556.
References
Note that CVR is reported as a percentage (Donahue et al., 2009; Lu et al., 2003; Taneja et al., 2020; Zhao et al., 2021) due to a dependency on the underlying baseline measures. That is, if baseline CBV () in vessel A is twice as high as in vessel B, there is double the vasoactive moiety in vessel A and double the blood volume change in absolute physiological units; thus, to ensure that the calculated CVR is the same in both vessels, values must be normalized to their respective baseline values.
A proportionality coefficient may be required if other units are chosen (e.g., when using mL for CBV, m/s for v, and mL/100 g/min for CBF).