The forward model is based on work by Uludag et al., 2009, originally developed to model gradient- (GRE) and spin-echo (SE) fMRI signals for magnetic field strengths between 1.5 and 16.4T (Uludag et al., 2009). As the fMRI signal change results from susceptibility variations, caused by dOHb in the case of neuronal activation or vascular challenges, the model can also be applied to DSC MRI utilizing any contrast agent that induces magnetic susceptibility changes. Here, this framework is generalized and applied to model 3T GRE signal changes induced by boluses of dOHb and Gd (please see the Discussion section for other DSC signal models and how they relate to the proposed model).

The signal model applied in this work was previously used to simulate fMRI signal changes in response to neuronally-driven blood oxygenation changes (Uludag et al., 2009). Given that the dOHb bolus used in this study is an externally driven susceptibility change of the same contrast agent as in fMRI, we applied the same biophysical model to DSC MRI. We generalized the framework to work for Gd contrast as well. To do so, the relationship between Gd concentration and susceptibility changes had to be determined. In doing so, we found that at 3T, a 1 mM increase of Gd is approximately equivalent in susceptibility to a 24.6% increase in dOHb (Eq. 8.2). Thus, Gd was modeled in the same way as dOHb, with differences only in the input bolus shape, duration, and susceptibility value (Supplementary Figures S1, S6). Consequently, any difference in calculated physiological parameters (e.g., CBV, CBF and MTT) can be attributed to differences in susceptibility and bolus properties.

Our model calculates the DSC MRI signal response to contrast agent in tissue, which results from both intravascular and extravascular relaxation rate changes in separate vascular compartments: artery, vein, capillary, arteriole, and venule. This has theoretical advantages over previous in vivo and simulation work which models DSC MRI through a “bulk model”, wherein there is no theoretical separation between extravascular and intravascular compartments (e.g., see Calamante et al., 2000; Patil and Johnson, 2011; Chappell et al., 2015). The contrast agent's influence on the proton spins in these compartments is not identical and should be treated as separate signal contributions (Zhao et al., 2007; Blockley et al., 2008; Uludag et al., 2009). The model described in Kiselev et al. incorporates extravascular and intravascular signal contributions and notes differences in the calculated and ground truth CBV and CBF of simulated voxels (Kiselev, 2001). A major difference between our work and Kiselev's is that the latter assumes an intravascular relaxation rate that is linearly dependent and an extravascular relaxation rate that is quadratically dependent on contrast concentration—our model assumes a linear extravascular and quadratic intravascular relationship.

Numerous experimental and theoretical papers report quadratic dependencies between contrast agent (dOHb or Gd) and the intravascular relaxation rate (Jensen and Chandra, 2000; Silvennoinen et al., 2003; van Osch et al., 2003; Zhao et al., 2007; Blockley et al., 2008). As a linear intravascular relaxation rate in our model would result in little to no dependency of the perfusion metrics on the susceptibility value, our experimental results are also a strong argument for a non-linear relationship between the contrast agent dose and intravascular relaxation rate. This is because the effects of contrast dosage, for example, would mostly cancel out in a linear relationship, resulting in little dependency. Regarding the extravascular relaxation rate, there is strong agreement based on rigorous Monte Carlo simulations that it is linearly related to the concentration of contrast agent (Boxerman et al., 1995; Kjølby et al., 2006; Uludag et al., 2009).

To the best of our knowledge, two studies conducted by Kjølby et al. (2006, 2009) are the only ones to account for linear extravascular and quadratic intravascular relaxation rates in DSC MRI simulations. The paper in 2006 (Kjølby et al., 2006) focused on deriving tissue relaxivity in response to Gd and did not quantify perfusion. The paper in 2009 (Kjølby et al., 2009) set out to quantify perfusion in response to changes in AIF partial volume. However, this work neglected separate vascular compartments in tissue. Therefore, even though components of our model have been proposed in some studies, the integration of all these components is novel, that is, incorporating quadratic intravascular and linear extravascular components from separate vascular compartments. In addition, the current work has investigated perfusion quantification more systematically and for a wider range of properties than previous studies.

Furthermore, it was not our aim to investigate perfusion quantification dependent on specific tissue properties, such as vessel branching, cellular packing, or vessel permeability; we instead focused on properties related to the applied contrast agents and intravascular volume fraction. Thus, we refer the reader to other DSC MRI signal models which have been developed to understand the effects of various tissue parameters (i.e., vessel heterogeneity, contrast agent extraversion) on signal and relaxation rate (Quarles et al., 2009; Semmineh et al., 2014; Digernes et al., 2017).

While Gd-induced signal changes are larger than those induced by dOHb, signal scaling is not homogeneous in all areas of the brain (Figure 4). When the ΔS_dOHb map is divided by the ΔS_Gd map, termed the ΔS ratio, we found a significant scaling difference between arterial, venous, GM, and WM regions. Specifically, large vessel voxels, on average, show a larger ΔS ratio than tissue voxels in both the simulations and experiments. The likely reason for this is due to non-linear signal scaling. When signal in the large vessels approaches zero, also known as saturation (Ellinger et al., 2000), signal change becomes non-linear with respect to dosage (doubling the dosage or CBV results in much smaller signal changes). Given that signal in response to Gd contrast can approach the noise floor, something that does not occur with a typical hypoxic bolus (Poublanc et al., 2021; Vu et al., 2021), there is an inherent difference in how signal is scaled between contrast agents with high vs. low susceptibility in large vessels. The scaling is more linear for dOHb and more non-linear for Gd, with respect to the amount of contrast agent.

Although to a much smaller degree, there is also a scaling difference between the GM and WM voxels, wherein the GM ΔS ratio between Gd and dOHb induced changes is lower than the WM ΔS ratio. This effect is observed in the experiments but not in the simulations (Figure 4). One potential reason for the observed difference between GM and WM might be macrovascular contamination (Chappell et al., 2015), particularly the so-called blooming effect (reviewed by Willats and Calamante, 2013). It is known that the administration of Gd results in signal contamination in the cortical GM, where magnetic field disturbances from large vessels extends to the surrounding tissue (reviewed by Willats and Calamante, 2013). Assuming that these vessels are large arteries, this effect would be largely scaled down when using dOHb contrast due to the far lower signal arising from the arterial vessels on the cortical surface (refer to Figures 2, 3). Further work would be required to isolate for and identify the magnitude of the blooming effect when using Gd and dOHb. In summary, the signal change induced by Gd is larger than that induced by dOHb, but the scaling is not uniform throughout the brain, confirmed both in the experiments and simulations. This is a novel finding, as no previous work has compared signal change maps obtained from contrast agents with vastly different molar susceptibilities.

Finally, the ΔS ratio in arterial voxels is lower than that in venous voxels for both the simulations and experimental data. The reason for this is more complicated than the prior findings and can best be understood by referring to the simulation findings in Figure 3A. Here, although results are displayed for ΔR2,max* as opposed to ΔS_max, the same rationale applies. The difference between Gd and dOHb is magnified in arterial voxels and is much less apparent in venous voxels—the division of ΔS_{max, dOHb} by ΔS_{max, Gd} therefore yields much lower values in arterial voxels. The mathematical reasoning for the difference between artery and vein will be discussed further in the following sections.

The simulations and experiments show a dependence of rCBV, rCBF, and MTT on user-controlled parameters, such as the type of contrast agent, choice of normalizing voxel, ΔS_aO₂, base-S_aO₂, and bolus duration. Figure 3 and Supplementary Figure 5 show that the area and magnitude of the arterial bolus relaxation curve at 3T do not increase linearly as the simulated blood volume increases, but rather increase and then decrease at higher simulated CBV occupancy. This counterintuitive pattern is also observed in previous work with Gd, where simulated arterial voxels with 80% CBV showed a greater signal change than those with 100% CBV (Kjølby et al., 2009). It is the predominant intravascular component of simulated arterial and venous voxels that results in a non-linearity in CBV dependence (Supplementary Figure 7); for simulated tissue voxels, this non-linearity is not observed due to the predominant linear extravascular relaxation rate. As such, the model predicts that normalizing tissue voxels to an arterial voxel with high blood volume results in higher calculated rCBV values than when normalized to an arterial voxel with an intermediate blood volume (Supplementary Figure 7).

As shown in Figures 5–8 and Supplementary Table 1, rCBV and rCBF attained with dOHb becomes more accurate and decreases as the bolus's ΔS_aO₂ increases and base-S_aO₂ decreases. Experimentally, this pattern should be more predominant when tissue is normalized to an arterial voxel with a high CBV, as from the theory and simulations, it is primarily the intravascular quadratic relaxation rate that is responsible for these dependencies (Supplementary Figure 7). The concept of the intravascular space contributing to perfusion overestimation is also acknowledged in previous works (Kjølby et al., 2009; Wirestam et al., 2010). As mentioned above, we show that increasing the susceptibility of contrast agent decreases the calculated rCBV and rCBF, without the confounds of saturation or partial volume errors in the case of our simulations. This pattern is more evident when using a lower susceptibility contrast agent, such as dOHb, which likely explains why previous researchers investigating the dosage of Gd found only a small effect of dosage on quantification (Alger et al., 2009; Wirestam et al., 2010).

Ultimately, the large deviations in the estimation of rCBV and rCBF between different dOHb and Gd paradigms are expected in accordance with the simulations. Of note, the rCBF values obtained in this work overestimate those from previous work using dOHb, which in turn overestimate physiological values determined by PET (Knutsson et al., 2010; Vu et al., 2021). This is similarly not unexpected as previous work applied a different experimental paradigm, which we now know from our simulations and validations will lead to a differing calculation of rCBF.

As an additional note, Bleeker et al. showed that it is beneficial to select an AIF adjacent to the intravascular space of an artery (Bleeker et al., 2009). Our work supports this claim for two reasons: firstly, our simulations show that normalizing tissue to an AIF with a low to intermediate simulated CBV results in higher accuracy of the calculated rCBV; secondly, our simulations reveal a lower dependence of the calculated rCBV on ΔS_aO₂ and base-S_aO₂ when tissue is normalized to a less intravascular AIF (Supplementary Figure 7).

From a clinical perspective, it might be beneficial to normalize tissue signal to a venous as opposed to arterial voxel, as the larger radius of the venous blood vessel results in less partial volume errors (reviewed by Willats and Calamante, 2013). Although there are still limitations to this type of analysis in DSC MRI (Knutsson et al., 2007), we sought to investigate how perfusion quantification depends on arterial vs. venous normalization.

In simulations and experiments, normalizing ΔR2,tissue*(t) to ΔR2,vein*(t) results in lower overestimations of rCBV, along with a reduced ΔS_aO₂ and base-S_aO₂ dependency in comparison to when ΔR2,tissue*(t) is normalized to ΔR2,artery*(t) (Figures 5, 6, and Supplementary Figure 5). As previously described, a lower base-S_aO₂ results in larger relaxation rate changes on account of the quadratic intravascular dependence on [Gd] and SO₂, which then influences perfusion estimation. Thus, if a vein and an artery are provided the same ΔS_aO₂ or [Gd] and are simulated with the same CBV, it is clear from the previous results that the voxel with a lower base-S_aO₂, as is the case in veins, will have a larger relaxation rate change (specifically when the voxel is mostly intravascular and when the contrast agent has a lower susceptibility or concentration). This theoretical discovery likely explains the high venous to arterial relaxation rate ratio observed in previous studies using dOHb in DSC MRI (Poublanc et al., 2021; Vu et al., 2021). When tissue signal is then normalized to the venous voxel with a larger relaxation rate change relative to artery, and thus, larger area under curve, a lower and more accurate estimate of rCBV is obtained. When normalizing to a vein, we found that varying the ΔS_aO₂ and base-S_aO₂ through different paradigms had less of an effect on perfusion quantification in both the simulations and experimental validations. A reduced contrast dosage dependency and lower calculated rCBV when normalizing to vein as opposed to artery is also observed in previous experimental results using Gd (Wirestam et al., 2010).

The experimental data reveal an overestimation of the calculated MTT in every participant when using dOHb in comparison to Gd for quantification (Figure 9). Vu et al. also reported high calculated MTT values when using dOHb as a contrast agent in a standard DSC MRI analysis (Vu et al., 2021).

The simulations reveal that bolus duration influences MTT and rCBF calculation (Figure 9). We also found that the thresholding of singular values causes a duration dependency—longer boluses resulted in the removal of more singular values from the AIF. It is known that as more singular values are removed, there is a greater underestimation of rCBF and, therefore, greater overestimation of MTT (Bjørnerud and Emblem, 2010). Thus, a Gd bolus with shorter duration and less removed singular values results in the calculation of a lower and more accurate MTT relative to the MTT calculated from a dOHb bolus. In our simulation work, when the dOHb bolus is shortened to the duration of Gd, the calculated MTT is more like that obtained with Gd (Figure 9). In previous work where Gd was doubled in concentration, an increase in the calculated MTT was found relative to a standard dose (Wirestam et al., 2010). Another study found that reducing the injection rate of Gd results in an increase to the calculated MTT (van Osch et al., 2003). As the bolus had increased in duration for these two studies, our findings support and help explain the observations.

Although noise was not simulated in our work, it is known that SNR is another major parameter influencing MTT calculation, wherein more singular values from the AIF matrix are removed when there is a lower SNR, resulting in a greater underestimation of rCBF and overestimation of MTT (Willats and Calamante, 2013). This may explain why the different dOHb paradigms in our work, although equal in duration, show an increasing estimation of MTT as the SNR decreases (Supplementary Table 1). Thus, a combination of low SNR and long bolus duration leads to higher estimations of MTT, supporting the finding that a dOHb bolus yields a higher MTT estimation than a Gd bolus.

For dOHb contrast to attain an SNR like that of Gd in a single bolus, the oxygen saturation must drop below physiological tolerance (Supplementary Figure 1 and Supplementary Table 1). However, DSC using dOHb has several advantages (Poublanc et al., 2021; Vu et al., 2021). Firstly, as dOHb is an endogenous contrast agent, it does not share the safety limitations associated with Gd (Schlaudecker and Bernheisel, 2009; Voth et al., 2011; Kanda et al., 2015; Lord et al., 2018; Rogowska et al., 2018). As well, [dOHb] does not saturate, recirculate, or accumulate, as the blood oxygenation is reset to a pre-determined value in the lungs using the RespirAct^®. Unlike Gd, the concentration of dOHb in the arterial blood can be determined from end-tidal gas analysis which may be useful as a precise AIF and aid in quantitative modeling (Ito et al., 2008; Fierstra et al., 2011). Finally, dOHb acquisitions can be repeated and averaged to improve the SNR (“Average dOHb” in Supplementary Table 1).

While the smaller hypoxic challenge (ΔS_aO₂ of 8%) administered in our study is equivalent to experiencing reduced S_aO₂ in moderately high altitudes (e.g., Aspen, Colorado), the larger challenge (ΔS_aO₂ of 24%) is closer to the S_aO₂ in the Everest Base Camp (Rojas-Camayo et al., 2018), albeit only for a short duration. Individuals may spend many months in these environments without additional oxygen or altitude sickness, minimizing the safety concern in our study where hypoxia only lasts for roughly 30 s.

We have not included T₁ effects in our simulations. According to Calamante et al., a TR of 1.5 s and flip angle of 73° would result in, roughly, a 10% underestimation of CBF relative to that calculated in the absence of T₁ effects, although this depends on the inflow rate (Calamante et al., 2007). This is to be considered in future studies but is not of major concern in our work as we are not investigating absolute values, but rather how relative quantification depends on various paradigms. It is not expected that the general patterns of the findings will change by including T₁ effects, especially in the absence of contrast agent extraversion. As well, it has been shown in previous work that there is no significant T₁ change in response to changes in dOHb, limiting the T₁ effect concern to Gd (Thulborn et al., 1982). A further limitation of our study and DSC MRI theory in general is the correct simulation of [Gd]_max. Studies have reported varying quantities for [Gd]_max (Patil et al., 2013; Kellner et al., 2018; Lind et al., 2020). Simulating the exact value of [Gd]_max is not a major concern in this work, but it should be noted that in vivo, a standard Gd dose likely results in a [Gd]_max that is larger than the simulated 2 mM used in the majority of our simulations.

In the model, we do not include contrast agent extravasation, which is known from previous modeling and experiments to elicit additional T2* and T₁ relaxation effects, limiting perfusion quantification accuracy (Quarles et al., 2005, 2009). Additional relaxation effects that result from leakage, dependent on the molecular size of contrast agent, can be implemented in future work that simulates pathological cerebral tissue. The limitation of neglecting this effect in our modeling is likely minimal as healthy subject tissue is assumed to have an intact blood-brain-barrier, with little to no leakage (Quarles et al., 2009).

Vasodilation from the hypoxic stimulus is another potential experimental limitation. Previous work showed that dropping S_aO₂ from 100 to 90% for a 20-min duration is the threshold for statistically significant cerebral vasodilation in humans (Gupta et al., 1997). According to more recent work, an S_aO₂ of ~70% results in increases to the CBF (Mardimae et al., 2012); however, the vascular response is delayed by about 3 minutes following hypoxia (Harris et al., 2013). We expect that there is no significant vasodilation in this work as the lowest S_aO₂ of 75% only lasts for approximately 30 s. Nevertheless, vasodilation would have to be accounted for in more drastic hypoxic paradigms, or in a diseased state, to avoid errors in quantification—this can easily be done with our model. Furthermore, for isocapnic hypoxia, as utilized in this study, no changes in the tissue oxygen metabolism or oxygen extraction fraction are expected (Ainslie et al., 2014; Gu and Jun, 2018).

As a future direction, dOHb contrast can be studied with a higher field strength scanner to further bolster SNR, which might be the best way to maximize the benefits of dOHb as long as the arterial signal does not become saturated. It is also known that at higher field strengths, the intravascular contribution in T2*-weighted acquisitions is heavily reduced relative to extravascular signal contribution (Uludag et al., 2009), which would lead to lower errors in perfusion quantification. It might be worthwhile, for a more detailed validation of the Gd to dOHb relationship determined in our modeling, to compare dOHb with low-dose Gd as contrast, where the low-dose of Gd better approximates the susceptibility induced by dOHb. We are currently in the planning stages of such a study. Another potential direction would be to investigate the use of dOHb contrast in clinical conditions, such as brain tumors, steno-occlusive vascular disease, and dementia. This could first be simulated with our framework and then validated with clinical scans. Finally, it would be useful to explore either a model-based DSC MRI analysis method or calibration approach using our DSC model, which should provide more consistent and accurate quantifications of perfusion.