Non-invasive perfusion MR imaging of the human brain via breath-holding

Dynamic susceptibility contrast (DSC) MRI plays a pivotal role in the accurate diagnosis and prognosis of several neurovascular diseases, but is limited by its reliance on gadolinium, an intravascularly injected chelated metal. Here, we determined the feasibility of measuring perfusion using a DSC analysis of breath-hold-induced gradient-echo-MRI signal changes. We acquired data at both 3 T and 7 T from ten healthy participants who engaged in eight consecutive breath-holds. By pairing a novel arterial input function strategy with a standard DSC MRI analysis, we measured the cerebral blood volume, flow, and transit delay, and found values to agree with those documented in the literature using gadolinium. We also observed voxel-wise agreement between breath-hold and arterial spin labeling measures of cerebral blood flow. Breath-holding resulted in significantly higher contrast-to-noise (6.2 at 3 T vs. 8.5 at 7 T) and gray matter-to-white matter contrast at higher field strength. Finally, using a simulation framework to assess the effect of dynamic vasodilation on perfusion estimation, we found global perfusion underestimation of 20–40%. For the first time, we have assessed the feasibility of and limitations associated with using breath-holds for perfusion estimation with DSC. We hope that the methods and results presented in this study will help pave the way toward contrast-free perfusion imaging, in both basic and clinical research.


Time course properties in bhDSC
The temporally averaged signal time courses were then converted to relaxation rate time courses (∆R 2 *(t)) by normalizing to the pre-and post-breath-hold baselines.
Figure 2 shows the breath-hold-induced ∆R 2 * time courses averaged across all subjects, for arterial (AIF), venous (VOF), gray matter (GM), and white matter (WM) voxels at 3 T and 7 T. ∆R 2 *(t) bolus dispersion and delay increases from input (AIF), to tissue (GM and WM), and ultimately to output (VOF), as is expected physiologically (Fig. 2).Unlike the other time courses, the AIF has a positive ∆R 2 * (to be discussed in Fig. 4 and the Discussion section).The bolus magnitude (defined here as the ∆R 2 *(t) integral (AUC )) was found to be ~ 1.5 times higher in GM relative to WM. Importantly, breath-holding at 7 T resulted in larger bolus magnitudes in the WM and GM (by ~ 2 times), artery (by ~ 2 times), and vein (by ~ 3 times) relative to those at 3 T-this is also illustrated in the subject-wise data (Figure S1) and simulation results (Figure S2).

Breath-hold CNR properties
The contrast-to-noise ratio (CNR) and associated GM-to-WM contrast were subsequently quantified (see Methods for details).Although eight breath-holds were averaged in the main results of our work, we also studied the effect of averaging fewer breath-holds to determine whether shorter scan time (fewer breath-holds) would yield sufficient CNR.
Figure 3 shows the subject-averaged CNR at 3 T and 7 T when modifying the number of breath-holds averaged (Fig. 3A,B), along with a subject-wise regression for CNR values between 3 and 7 T (Fig. 3C).CNR and GMto-WM contrast were significantly higher at 7 T relative to 3 T (p CNR = 0.0012; p GM-to-WM = 0.0009) (Table 1).As more breath-holds were averaged, CNR increased significantly for both 3 T (p = 0.00006) and 7 T (p = 0.00007); at both field strengths, we found that CNR gains fit to a radical ( a • √ b • x + c ; R 2 = 0.99) as a function of boluses averaged (Fig. 3B).Of note, the subject-wise regression of CNR values between field strengths (Fig. 3C), which provides some insight into subject repeatability, was very high (R 2 = 0.7).According to the regression slope, CNR values were 1.45 times higher at 7 T relative to 3 T.

Novel determination of an AIF for bhDSC
A prerequisite for DSC MRI is the presence of signal change upstream of tissue (typically at the arterial level), which is then used as an input function (i.e., AIF) for perfusion quantification in the tissue and veins.In Gd-or hypoxia-based methods, arterial signal change results from an increase in paramagnetic contrast agent in the artery, which subsequently passes through to the tissue 6,8 .However, given that blood oxygenation in healthy subjects is typically fully saturated in the major cerebral arteries, an increase in blood flow induced by CO 2 following breath-holding is not expected to significantly change the amount of dOHb in the arteries [26][27][28] .Consequently, arterial magnetic susceptibility is not expected to change much as a result of breath-holding.However, hypercapnia uniquely results in vasodilation 29 -therefore, while MRI signal is not expected to change in the arteries from paramagnetic contrast agent, we expect that arterial vasodilation will result in a signal decrease, particularly at higher magnetic field strength (see Discussion).Thus, we developed a novel, alternative framework for selecting the AIF when using breath-holds (or hypercapnic paradigms in general).
Figure 4 shows voxels exhibiting a negative signal change (i.e., a positive relaxation change), and, following averaging of these voxels in the middle (MCA), posterior (PCA), and anterior cerebral arteries (ACA) (refer to Methods for selection criteria), displays the subject-averaged AIF time courses at 3 T and 7 T in response to breath-holding.Voxels with a positive relaxation rate change are located adjacent to the ventricles and in regions containing/adjacent to veins or arteries-the magnitude of the associated relaxation rate changes are significantly higher (p < 0.0007) at 7 T in comparison to 3 T.The presence of voxels displaying a positive relaxation rate change is, at first glance, counterintuitive given that the predominant effect of breath-holding is a decrease in dOHb (i.e., decrease in the relaxation rate).However, when we simulated (Figure S2) a voxel containing a large vessel, with little-to-no change in dOHb in combination with vasodilation, a substantial relaxation rate increase was in fact observed, in agreement with the experimental results.Thus, the simulations and experimental findings support the hypothesis that arterial vasodilation is responsible for the substantial relaxation rate increase observed in the arteries as a result of breath-holding.

bhDSC perfusion measurement
The AIF was first scaled by the integral of the VOF, composed of voxels in the superior sagittal sinus (SSS)-this step was conducted as, like in tissue, the VOF integral is representative of the amount of contrast agent (dOHb), whereas, unlike in tissue, the AIF integral is representative of the degree of vasodilation.Following this scaling step, CBV, CBF, and MTT maps were calculated using a standard, truncated singular value decomposition (SVD) analysis, with an SVD noise threshold of 20% and hematocrit correction factor of 1.45 11,[30][31][32][33] .Please note that the reported perfusion values are considered relative, as these values are scaled by some unknown set of factors which depend on acquisition, analysis, and tissue parameters (see Methods) 11 ; therefore, we report the CBV and CBF values as (a.u.).As will be discussed, the scaling is not random but fully determined by acquisition and analysis parameters.
Figure 5 shows the subject-averaged CBV and CBF perfusion maps (Fig. 5A), subject-wise GM and WM perfusion values (Fig. 5B), and voxel-wise regressions for CBV and CBF between 3 and 7 T (Fig. 5C).At 3 T, CBV values are slightly higher (p = 0.108) and CBF values are significantly higher (p = 0.0027) than at 7 T (Fig. 5B).The attained GM-to-WM perfusion contrast (i.e., ratio of GM-to-WM perfusion values) was significantly higher at 7 T relative to 3 T for both CBF (p = 0.018) and CBV (p = 0.016) values.In addition, as shown in the subject-wise CBF maps (Figure S3), clear GM-to-WM contrast is observed in all subjects at 7 T, but not at 3 T. High coefficients of determination are associated with the voxel-wise linear regression of 3 T and 7 T CBV (voxel-wise  www.nature.com/scientificreports/R 2 = 0.53) and CBF (voxel-wise R 2 = 0.58) values (Fig. 5C).This regional agreement is further supported by GMnormalized CBV and CBF maps (Figure S4); more pronounced differences can be observed in the WM, although these values are not reliable in bhDSC or ASL due to lower CNR and longer vascular transit times.Of note, MTT values in GM (Figure S5) are slightly higher, although not significantly (p = 0.74), at 7 T as opposed to 3 T (MTT 7T,GM = 7.36 ± 2.52 vs MTT 3T,GM = 6.96 ± 2.11).Due to low CNR, MTT values are not reported from the WM.

Comparing bhDSC with ASL
For validation, we compared our bhDSC results with 3 T ASL data obtained from the same subjects (refer to Methods for calculations of CNR and CBF for ASL).Note that CBF values for ASL are not absolute, as a calibration scan was not obtained-thus, only a relative comparison is conducted for CBF.

Effect of vasodilation on CBV estimation
Given that bhDSC relies on dynamic vasodilation for bolus generation in the tissue, unlike in other methods such as ASL or Gd-DSC, it is unclear as to whether baseline perfusion measures can truly be estimated during a vasodilatory challenge.To address this concern, we developed simulations (refer to Supplementary Materials) 11 to assess the magnitude by which tissue vasodilation introduces error into calculated baseline perfusion measurements.
Figure 7 shows the effect that simulated tissue vasodilation has on CBV quantification.According to the simulations, the induced vasodilation that occurs during a breath-hold results in CBV underestimation, relatively independent of baseline CBV values and tissue composition.Given an approximate relative vasodilation (ΔCBV) of 4-9% during breath-holding, the associated CBV underestimation is ~ 20-40%.Thus, for example, a CBV of 4% would be underestimated to be around 2.8%.Note that a range is provided as, depending on breath-hold performance, the amount of vasodilation will vary between subjects and, consequently, the underestimation of CBV values will be subject dependent.Although globally underestimated, it appears that regardless of tissue blood www.nature.com/scientificreports/volume or vessel composition, underestimation does not vary much and, therefore, regional scaling differences are not expected.In other words, even though the absolute quantitative values are affected by vasodilation, the relative distribution of values remains largely unaffected.However, in the event that different voxels have substantial variations in vasodilatory capacity (reflected by more variability in the x-axis for a given subject), it can be expected that regional scaling differences due to vasodilation will be observed (see Discussion).Although not shown here, CBF underestimation from vasodilation shows the same behavior as CBV underestimation; thus, the MTT is not notably affected.

Discussion
For the first time, we performed a DSC-MRI perfusion analysis using a breath-hold task, without exogenous contrast or additional medical device equipment.We have found the following: • The AIF during a breath-hold task is uniquely characterized by a negative signal change, likely representing vasodilation, as supported by both simulations and the literature.• bhDSC-calculated CBV and CBF are generally within the physiological range of values reported in the lit- erature using established DSC-MRI approaches.• bhDSC-calculated CBF demonstrates regional congruency and voxel-wise linear agreement with CBF deter- mined using ASL.• At 7 T, the breath-hold task yielded significantly higher CNR (p = 0.0012) and GM-to-WM contrast (p = 0.0009) relative to 3 T and ASL (p = 0.0043 and p < 0.0001, respectively).• Tissue vasodilation, unique to bhDSC, yields a global CBV and CBF underestimation of ~ 20-40%, but assum- ing relatively homogeneous vasodilatory capacity, error is consistent across voxels with different vascular properties.

From breath-hold to signal change
It is well established that hypercapnia induced by breath-holding produces a quantifiable T 2 * signal increase in the tissue, which has been particularly useful in CVR and calibrated BOLD studies 24,25,34,35 .Of note, mild hypoxia is also known to accompany hypercapnia during breath-holding, although, given a ~ 0.65% reduction in arterial oxygen saturation during a 16 s breath-hold 27,36 , the impact of breath-hold-induced hypoxia on T 2 * signal change is expected to be negligible in comparison to the effect attributed to hypercapnia.It has long been known that increased P a CO 2 (arterial partial pressure of CO 2 ) leads to increased blood flow in the arteries, and subsequently, in tissue capillaries and veins 37 .This occurs either directly through CO 2 acting on the smooth muscle cells surrounding the arteries/arterioles, or indirectly through CO 2 acting on the vascular endothelium, in both cases resulting in vasodilation and increased blood flow 21 .Therefore, a breath-hold, which yields a transient rise and fall (i.e., bolus) of P a CO 2 , is a relatively simple way to induce a transient bolus of CBF.
Of note, the breath-hold duration has a sigmoidal relationship with the resulting change in P a CO 2 36 , which then also has a sigmoidal relationship with the resulting change in CBF 38 .Given that oxygen metabolism is assumed to be constant during mild hypercapnia 23,26,28,39,40 , the CBF bolus resulting from arterial/arteriolar/capillary vasodilation now permits the cerebral tissue to be supplied with more oxygen than it requires, leading to an increase in blood oxygenation and a reduction in dOHb, in accordance with the bolus shape and magnitude of P a CO 2 in tissue.
We refer to what has just been described as the CVR hypothesis, wherein the observed arterial, tissue, and venous signal changes are a result of CO 2 -induced vasodilation throughout the vascular tree.That is, the increased CO 2 following breath-holding leads to an increase in arterial blood volume, in turn resulting in increased flow.In this hypothesis, it is predominantly CO 2 -induced vasodilation/flow increases prior to the tissue which result in an oxygenation bolus at the tissue level.However, a second mechanism-the cardiac hypothesis-may alternatively describe the observed oxygenation bolus.Here, the previously described CVR effect may occur in addition to a cardiac output increase, which propagates throughout the cerebral vasculature as a bolus of increased blood flow (in accordance with the breath-hold duration).According to the literature, the likelihood of a cardiac hypothesis is very low given that cardiac output can increase, decrease, or remain constant depending on the breath-hold maneuver 41 .In fact, Sakuma et al. found that large lung volume breath-holding (end-inspiration) results in a reduction of cardiac output, which would theoretically yield the reverse effect of that observed in our data 42 .
Nevertheless, these mechanisms predict a change in paramagnetic dOHb concentration in the tissue vasculature and draining veins, which is the basis for T 2 * signal change in fMRI [43][44][45][46] .The magnitude of these signal changes scale with the underlying baseline CBV, allowing for measurement within a DSC framework.In the current study, we observed a ~ 10-15 s delay between the initiation of hypercapnia and the resulting T 2 * signal change bolus onset (Fig. 1), which is similarly observed in previous hypercapnia studies and is attributed to the time for CO 2 build-up in the lungs and the lung-to-brain travel time of blood 47,48 .Ultimately, the decrease in dOHb results in a transient increase in MRI signal, and the magnitude of the signal change is dependent on the magnitude of ΔP a CO 2 , field strength, pulse sequence, and tissue composition 49 .

Determination of an AIF for bhDSC
To measure CBV, CBF, and MTT using DSC MRI, the tissue bolus must be deconvolved with an input bolus, ideally from the cerebral arteries (AIF), that acquires dispersion and delay as it travels to and through the tissue, in accordance with indicator dilution theory 6,8 .While an AIF is expected and routinely determined when using Gd or hypoxia contrast, the presence of an AIF and its associated properties are unknown in bhDSC.To address this, we examined the physiological processes underlying arterial signal change as a result of breath-holding.As argued above, dOHb levels in healthy subjects do not change much in the arteries during breath-holding as there is effectively no oxygen exchange at this level of the vasculature and arterial blood oxygenation is already near-saturation 26,28 .Thus, given that most voxels in a typical T 2 *-weighted MRI acquisition (2-4 mm voxel dimensions) likely consist of more than one tissue/vessel type, any observed signal increase in and around the MCA, for example, may represent signal changes from surrounding, early-perfused cortical tissue and nearby veins, and using this as an AIF would violate basic tenets of indicator dilution theory.
Experimentally, we noticed that the MRI signal, which increased throughout most of the brain in response to a breath-hold, decreased in many voxels containing large arteries and/or veins, particularly when adjacent to CSF, with a larger observed effect at 7 T (Fig. 4).To better understand this phenomenon, we modeled (Figure S2) an arterial voxel with vasodilation (i.e., CBV change) but no dOHb change.In doing so, we found that vasodilation resulted in a sizable signal decrease (i.e., apparent relaxation rate increase) in accordance with our experimental findings.The reason for this phenomenon relies on the fact that the intravascular transverse relaxation rate is substantially higher than the extravascular tissue and CSF relaxation rates at 7 T, and to a much lesser extent, at 3 T, for which the intravascular transverse relaxation rate is higher than CSF relaxation 50 .Thus, an increase in CBV with little-to-no change in blood oxygenation leads to an MRI signal decrease (Fig. 4C).This signal origin has also been hypothesized as an explanation for the initial signal dip observed in fMRI studies 49 and has been posited in other work 51,52 .Although its magnitude depends on vessel orientation and voxel composition, measurable signal change is expected from our simulations to occur in any voxel where a relatively large vessel increases in volume with negligible change in blood oxygenation, displacing extravascular tissue and/or CSF, particularly at 7 T. Of note, this means that tissue voxels co-localized with larger vessels and CSF may display both an initial signal decrease from vasodilation and a subsequent signal increase from the resulting dOHb change.
Our findings are also supported by the literature, where negative signal change during hypercapnia was found in and adjacent to the ventricles 47,53 .At 3 T, this effect results almost entirely from the difference between CSF and intravascular signal, as previously described 53 .Voxels characterized by a signal decrease were observed in higher abundance at the arterial and venous level in a previous fMRI study at 7 T 54 , as is similarly observed here.www.nature.com/scientificreports/While we were unable to sample P a CO 2 in the cerebral arteries while directly imaging changes in vascular occupancy-a potential avenue for future work-there is very strong experimental and theoretical evidence that the apparent positive T 2 * relaxation rate change at the arterial level represents a time course of arterial vasodilation in response to changes in P a CO 2 .This phenomenon is not expected in Gd or hypoxia studies, where the relaxation rate is known to increase in the artery due to increased paramagnetic contrast agent, and not vasodilation.Although it is not the traditional input of contrast agent-induced signal change, arterial vasodilation does in fact represent the input induced by hypercapnia, likely reflecting the input P a CO 2 time course in the brain, which then passes through and acts on downstream vasculature to yield tissue signal change in accordance with the shape and magnitude of the P a CO 2 bolus.The careful selection and averaging of voxels in and around the MCA, PCA, and ACA exhibiting this apparent positive relaxation rate change provided us with a novel solution for deconvolving the tissue data with an input function.The fact that arterial signal change is opposite in sign to signal change observed in much of the brain will likely make it easier to identify an AIF using bhDSC in comparison to traditional DSC techniques.Again, given that we have not acquired P a CO 2 time courses in our data, future work would need to establish the quantitative relationship between P a CO 2 and MRI signal change in and around the arteries.
Of note, we recognize that vasodilation, and thus negative signal change, may also arise at the arterial level from other sources, such as changes elicited in the cardiac cycle.However, we are confident that these extraneous effects are either aliased due to the lower temporal resolution (TR = 2 s) or averaged out during the process of breath-hold bolus averaging.

CNR in bhDSC
There is a significant increase in CNR and GM-to-WM contrast as field strength increases (Fig. 3; Table 1).These findings are supported by our simulations (Figure S2) and studies that have demonstrated a higher transverse relaxivity of dOHb as a function of field strength 49,55,56 .Thus, for the same breath-hold duration and performance, higher signal change will be observed at 7 T relative to 3 T.In fact, 1 min of scan time at 7 T resulted in approximately the same CNR and GM-to-WM contrast as ~ 4 min of scan time at 3 T.This supports the notion that a bhDSC analysis will yield more robust and, hence, clinically valuable results at higher field strength, within a shorter scan duration.Although there is also a significant increase in CNR as the number of breath-holds averaged increases, we found the smallest CNR and GM-to-WM contrast gains when comparing 6 versus 8 breath-holds averaged (Fig. 3A,B), in agreement with the known CNR ∝ √ Averages relationship.Our data sug- gest that bhDSC (TA: 8.5 min) yields superior CNR with respect to ASL (TA: 6.5 min), particularly at 7 T where CNR is ~ 1.6 times higher than ASL CNR.Of note, bhDSC CNR at 7 T is even significantly higher than ASL CNR (p < 0.05) when using 6 breath-holds (TA: ~ 6.5 min).As both methods are non-invasive, this is an important finding given that CNR is a known limitation of ASL, and bhDSC can yield additional perfusion parameters (CBV and MTT) in comparison to standard ASL.
The minor increase in WM CNR as a function of field strength requires closer examination.This finding may be partially explained by the fact that fMRI signal in large pial veins is known to increase more than linearly as a function of field strength (and these veins are often co-localized to the GM mask).Thus, we expect that to some extent, GM CNR will increase more than WM CNR as a function of field strength, leading to a difference in the GM-to-WM ratio.It is also possible that the calculated GM-to-WM perfusion ratio varies as a function of field strength for reasons other than this partial volume effect, although this is unclear.Until validated further, bhDSC may be most useful for studying GM perfusion and large vessel flow.

bhDSC perfusion measurement
The subject-averaged CBV and CBF values and maps (Figs. 4 and S3; Table 1) are generally within the documented range of previously reported values (CBV GM = 3-11 mL/100 g; CBF GM = 52-137 mL/100 g/min) 30,57,58 .While high GM-to-WM contrast is observed in all individual subject maps at 7 T, this is not always the case at 3 T (Figure S3), suggesting that bhDSC will yield perfusion maps that are more reliable at higher field strength.
The discrepancy between 3 and 7 T values supports the notion that DSC measurement is relative, with a dependence on field strength, amongst many other parameters.The high-level explanation for the CBV and CBF measurement discrepancy is that tissue AUC increases less than VOF AUC as field strength increases (Figs. 2 and  S2).Specifically, we found that AUC GM nearly doubles while the average AUC VOF nearly triples as field strength increases from 3 to 7 T. Given that CBV is calculated as a ratio of AUC GM to AUC VOF (and CBF is scaled by this ratio), the aforementioned scaling difference results in lower calculated values of CBV and CBF at higher field strength-this is also supported by our simulations (Figure S2).This scaling difference likely results from the non-linear intravascular contribution, which is a larger signal component at lower field strength (approximately 10-40% of the signal at 3 T in comparison to about 0-5% of the signal at 7 T in accordance with Uludag et al. 2009) 49 .To reiterate, as there is no commonly accepted pipeline for absolute DSC MRI quantification due to the many quantification dependencies that exist, calculated perfusion values are, to-date, considered relative, and only these relative values are considered in clinical care.Nevertheless, there is a strong voxel-wise linear correlation between 3 and 7 T tissue perfusion values (Fig. 5C) and visual agreement in the GM (Figure S4), indicating that there are no major regional gray matter scaling differences between field strengths in the GM of healthy subjects.We hope that our documentation of the perfusion estimation differences between 3 and 7 T will aid in the development of corrective measures in the future, for absolute quantification.
As previously described, the CBF maps are generally well correlated between bhDSC and ASL, with the largest deviations observed in surrounding vascular territories (Fig. 6C).The observed incongruence between ASL and bhDSC is expected in the arterial and venous territories surrounding the GM: methods such as bhDSC and Gd-DSC track a contrast agent (dOHb and Gd, respectively), which does not traverse the blood brain Vol:.(1234567890 www.nature.com/scientificreports/barrier in healthy subjects, yielding large signal changes in the capillaries/tissue, supplying arteries, and draining venous vasculature.In methods such as ASL and H 2 O-PET, the tracer is labeled water, which in humans, largely exchanges with the tissue, and therefore only yields high signal change in the capillaries and tissue (low signal change in the arteries/veins).Thus, regional incongruency in territories containing large vessels is expected between bhDSC and ASL.Although we did not acquire data from multiple sessions per subject at a given field strength, we found that when bootstrapping different bolus combinations to calculate GM CBV at 7 T (i.e., averaging 4 boluses for a total of 70 bolus combinations), the intraclass correlation (ICC = 0.63) was in the moderate to substantial range (Figure S7), with the lowest precision observed in subjects where the breath-hold bolus shape was not well-defined (subjects 9 and 10).In addition, we observed a strong correlation in subject-wise CNR between field strengths (Fig. 3C).Both of these findings support the notion that bhDSC, at the very least, is a moderately repeatable method.
While MTT GM values (Figure S5; Table 1) are twice as high as values obtained in studies using Gd contrast, they are similar to, if not lower than those reported in studies using hypoxia contrast 11,15,30,58 .Once again, this discrepancy likely results from another known relative DSC quantification dependency, namely: we have previously found longer bolus durations to result in higher calculated MTT values when using a standard DSC MRI framework 11 .A breath-hold of 16 s induces, on average, a 35 s signal bolus, which is substantially longer than a Gd bolus (but similar in duration to a typical hypoxia bolus), and consequently yields a higher calculated MTT due to the inherent limitations of singular value decomposition.Although slightly higher at 7 T, MTT measurements are not statistically different between field strengths.Of note, in accordance with the central volume principle 31 , the higher calculated MTT values explain the slight underestimation of CBF in comparison to average values reported in the literature 30,57,58 .
Finally, the delay maps are useful for addressing whether the breath-hold induces a bolus that traverses the brain in a manner dependent on the vascular hierarchy/topography, acquiring delay/dispersion as it progresses.According to our delay maps (Figure S8), the bhDSC method yields delay values across the brain in accordance with the ordering of bolus arrival times reported in the literature when using Gd as a contrast agent 59 and breathholding for CVR analysis 60 .This provides additional evidence that bhDSC involves a bolus whose delay and dispersion properties are dependent on the location within the cerebral vascular network where tissue is perfused.

Clinical and physiological considerations for bhDSC
There are notable benefits associated with bhDSC in comparison to traditional DSC methods, including the lack of exogenous contrast and/or medical device equipment, and their associated limitations (i.e., patient discomfort, contrast agent extravasation, medical device maintenance, and additional costs) 13,61 .Given that multiple perfusion parameters can be acquired without the necessity of (potentially) scanner-restrictive multi-delay acquisitions, there are also notable benefits with respect to ASL.However, it is important to consider whether there are physiological and/or physical processes, unique to breath-holding, which limit bhDSC's clinical utility in comparison to available state-of-the-art methods.
In a standard DSC analysis 32 , it is assumed that CBV and CBF do not change, a condition which is fulfilled during mild hypoxia or Gd.In bhDSC, signal change induced by a dynamic vasodilation time course in tissue may confound the calculated perfusion metrics.To investigate, we simulated (see Supplementary Materials for simulation framework) a 16 s breath-hold in addition to vasodilation in tissue and found global CBV and CBF underestimation by 20-40% (Fig. 7).We found the underestimation to be almost identical in magnitude for voxels with differing vascular compositions and CBV values, thus, relative values are expected to be preserved across the brain.On the other hand, given the presence of other global quantification errors in DSC, global scaling underestimation from vasodilation is not expected to hinder the utility of this method, particularly if a global correction factor is introduced.Thus, from our simulations, vasodilation in the tissue is not expected to confound the relative distribution of CBV values in healthy subjects.However, this may not be the case in patients, for which there are more regional variations in vasodilatory capacity compared to healthy subjects.Here, it may be useful to develop and employ correction strategies for the relative perfusion values in diseased tissue when using bhDSC.
Another consideration is the potential variability in the oxygen metabolism rate (CMRO 2 ) during breathholding, as this would additionally violate an important DSC analysis assumption.There is very compelling evidence that shows a lack of significant CMRO 2 variability during mild hypercapnia or breath-holding, especially for the short, 16 s breath-hold duration employed in our study 23,26,28,39,40 .
Causality of flow in the vascular network is another important consideration.In healthy subjects, the cerebral vasculature is assumed to respond relatively homogeneously to a vasodilatory stimulus 38,62 .However, in certain patients (i.e., those with steno-occlusive disease and/or vascular steal physiology), the affected arteries may have reduced vasodilatory capacity, which would result in tissue downstream from the affected vessels yielding reduced calculated CBV values in comparison to tissue supplied by other vessels (and potentially negative values in the case of vascular steal) 63 , even though the ground truth blood volume values may not be reduced.In these patients, apparent hypoperfusion relative to the rest of the brain may be a result of combined hypoperfusion and reduced vasodilatory capacity in the supplying vasculature, limiting bhDSC's validity in this cohort without any additional data or correction algorithms.
Finally, unlike Gd-and hypoxia-based methods, the signal change in bhDSC reflects an oxygenation bolus that only begins at the tissue level-thus, signal change in bhDSC is capillary/venous dominated.When comparing bhDSC with Gd-DSC perfusion results, it is expected that there will be some incongruency in brain regions which have substantially more arterial blood.However, this difference also creates an opportunity to combine information from multiple methods to isolate different blood pools.www.nature.com/scientificreports/

Limitations and future directions
Given that dOHb has a far lower molar susceptibility than Gd 11 and ΔdOHb is relatively small during breathholding 26,36 , it can be expected that bhDSC will yield substantially lower CNR than Gd-DSC.However, as illustrated in this work, both the accuracy and precision of perfusion measurement with bhDSC is still sufficient on an individual subject level at 7 T and on a group level at 3 T.At lower field strengths (such as 3 T), incorporating physiological and physical noise correction methods and acquiring multi-echo T 2 * acquisitions 64 , may enable bhDSC to yield more robust perfusion metrics for individual subjects.The breath-hold itself confers some limitations in comparison to ASL.Patient compliance is required during the repeated breath-holds of 16 s, which may be difficult for certain subjects, particularly the elderly, and those with neurological or lung disease 65,66 , However, breath-holds are often applied in patients for body imaging, and as previously described, sufficient contrast can be obtained at 7 T with fewer breath-holds, thereby, reducing the duration required for patient compliance.Another limitation is that movement associated with breath-holding may result in detectable motion during the scan.However, in the present study, we found this to be minimal in effect and correctable using post-processing (motion correction).As well, the breath-hold paradigm requires either a projector setup in the scan room (and the associated programming), or auditory cues through headphones.Finally, there is inter-subject variability in the breath-holds, which is reflected in the variability of bolus shape and magnitude between subjects (Figure S1).The variability (i.e., coefficient of variation (CV)) of our calculated perfusion metrics (Table 1) is on par with that observed in a previous study 15 which used a gas control system for DSC perfusion measurement at 3 T (CBV CV = 0.33; CBF CV = 0.42), limiting the concern that the contrast bolus needs to be standardized across subjects to ensure precision of the calculated perfusion values.Nevertheless, precision should be enhanced in future work by optimizing the breath-hold approach-providing more detailed instructions to the subjects, including whether an end-inspiration or end-expiration is employed, and potentially utilizing longer rest durations between breath-holds 67 .
While we have performed a direct comparison with ASL, a direct subject-wise bhDSC vs Gd-DSC comparison will allow for a stronger validation of the technique's utility and limitations in perfusion imaging.As well, the clinical utility of the bhDSC technique should be validated in different patient cohorts, to determine whether all or a subset of those with vascular brain pathology will benefit from bhDSC alone or in combination with additional perfusion imaging techniques (i.e., ASL).Finally, to provide more clinical centers with access to this method, it will be prudent to assess the viability of and optimize the bhDSC technique at lower field strengths.

Conclusion
For the first time, we leveraged breath-holds to measure perfusion with a DSC analysis.This was made possible by developing and applying a novel arterial input function strategy, tailored for breath-hold and hypercapnic paradigms.In doing so, we found that perfusion values (i.e., CBF and CBV) yielded higher CNR and GM-to-WM contrast at 7 T, were within the range of physiological and literature-reported values, and demonstrated high regional correlation with ASL CBF data obtained on the same subjects.While bhDSC will need to be further validated with other perfusion techniques in healthy subjects and in various diseases, we are convinced that our findings will aid the implementation of contrast-free perfusion imaging, in both basic and clinical research.

Subjects
This study 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 (age: 30.4 ± 9.4 years, 3 female).Please note that the same subjects were scanned at both 3 T and 7 T.

MRI sequences and experimental protocols
MRI data were acquired at the Center for Neuroscience Imaging Research at Sungkyunkwan University on the Siemens 3 T Prisma (Siemens Healthineers, Erlangen, Germany) and Siemens 7 T Terra (Siemens Healthineers, Erlangen, Germany) using the commercially available 64 channel head/neck and 32 channel head coils (Nova Medical, Wilmington, USA), respectively.

Structural MRI
The 3 T parameters were as follows: 3D-MP2RAGE 68  ASL MRI ASL scans were only acquired at 3 T.The ASL data were acquired with 2.0 mm isotropic spatial resolution, FAIR QUIPSS II labeling scheme and 3D gradient-and-spin-echo (GRASE) readout 69  measurements of opposite phase-encoded (A-P) data were also acquired for distortion correction.Please note that the spatial resolution has not been optimized for 7 T in this study.We have instead chosen the same spatial resolution as at 3 T to allow for an easier and more direct comparison.
For the breath-hold scans, subjects were instructed to fixate on a black cross at the center of an iso-luminant gray screen and a countdown timer, which indicated when to breathe regularly and when to perform a breathhold.Specifically, each breath-hold block was composed of three sections: 10 s of preparation with regular breathing, 16 s of breath-holding, and 34 s of regular breathing (Fig. 1).This was repeated nine times to obtain multiple boluses for averaging, with an added final baseline of 80 s.Only the first eight boluses were used for each subject as the response to the final breath-hold bolus was not fully included during the acquisition for some subjects.
The MP2RAGE structural data were pre-processed using presurfer (https:// github.com/ srika sh/ presu rfer) and skull-stripped using SynthStrip 72 .The background denoised T1-weighted UNI image was segmented using FSL FAST 73 to obtain three tissue classes corresponding to gray matter (GM), white matter (WM), and cerebrospinal fluid (CSF).The breath-hold data underwent slice-timing (3dTshift, AFNI), motion (mcflirt, FSL), and distortion correction (topup, FSL) 74,75 .The pre-processed breath-hold data were linearly detrended and temporally filtered by averaging each signal time point with a 1 × 5 Gaussian kernel.The first eight breath-hold blocks were then averaged, wherein each of the eight "boluses" was allotted a 72 s window centered at the bolus maximum and averaged-an example of an original and averaged time course is shown in Fig. 1.The percentage maximal signal change (∆S) and ∆S-to-noise ratio, hereafter termed contrast-to-noise ratio (CNR), were then calculated from the bolus-averaged signal time course: In Eq. 1, S 0 is the average signal of ten temporal volumes before and after (defined as the time that the signal returns to the average pre-bolus baseline) the bolus and S Max is the maximum signal increase.The CNR is computed as shown in Eq. 2, where ε t is the standard deviation of signal for ten total temporal volumes before and after the selected bolus.GM-to-WM contrast (GW contrast ) was also calculated as follows: where CNR GM and CNR WM represent the mean GM and WM CNRs, respectively.S(t) was then converted to the change in transverse relaxation rate time course ( �R * 2 (t)): The above equations were applied to T 2 *-weighted breath-hold data from all subjects at 3 T and 7 T to generate S , CNR , and �R * 2 (t) maps.Note that changes in intravascular or extravascular volume/occupancy, during vasodilation for example, will result in a change in the 'apparent' R * 2 (t).The ASL data first underwent motion correction (mcflirt, FSL), resampling to anatomical space (nearest neighbour), and masking (using the skull-stripped anatomical mask) 75 .Pairwise subtraction across the four tag-control pairs and subsequent averaging was performed to calculate CBF maps.Given the lack of a calibration scan, the CBF maps are not absolute, and only a comparison of relative CBF maps was conducted between ASL and bhDSC.ASL CNR maps were also calculated based on Eq. 2, where S Max is the average label signal, S 0 is the average control signal, and ε t is the standard deviation across the controls.
To facilitate group analyses, these maps, along with the later described bhDSC perfusion maps, were transformed to FSL's MNI152 2 mm space.The processed bhDSC data were first registered to the subject's anatomical space (6 dof, FSL flirt) and the anatomical data were then non-linearly registered (fnirt, FSL) to the FSL MNI152 2 mm template 71,76 .The two transformation matrices were then combined into a subject-specific native-to-MNI warp and applied to all subject-specific ASL and T 2 * maps.The MNI space-transformed maps were then averaged

Gray and white matter segmentation
GM and WM masks were generated, at both 3 T and 7 T, using whole-brain anatomical data (FSL FAST).The partial volume estimate maps were then thresholded at 0.9, binarized, and transformed from anatomical to GRE-EPI space for each subject (nearest-neighbor interpolation).Arterial and venous voxels were removed from the masks by thresholding the GRE-EPI temporal standard deviation map (only voxels with values in the lowest 10% of the ε t range were kept in the mask) and using the output to mask the binarized FAST segmentations-the assumption being that voxels with high CBV have higher physiological noise due to pulsatility, leading to higher values of ε t 11,77 .The whole-brain GM and WM masks were then used to obtain average GM and WM values (Table 1) from each subject's ASL, bhDSC perfusion, ∆S, and CNR maps.GM and WM masks were also generated in standard space from the FSL MNI152 T 1 2 mm template using the procedure described above.

Defining the AIF and VOF
In DSC MRI, it is critical to define an input function that is then used to normalize and deconvolve the tissue data-this input is defined at the level of the artery and is consequently termed the arterial input function (AIF) 61 .Typically, to represent the AIF, arterial voxels are identified with S values above the 95 th percentile and zero delay (see delay measurement below).However, given the spatial resolution of GRE-EPI (typically 2-4 mm in each dimension), signal change derived at the arterial level may be contaminated by signal changes from cortical tissue and nearby veins.With this in mind, we developed an alternative, novel method to obtain an AIF from the breath-hold data, which can be used in any study that uses hypercapnia to calculate perfusion at high or ultra-high field strength.In essence, the selection criteria are the same as those described above, but voxels yielding negative S values, with magnitude above the 99th percentile (S Max is replaced with S Min in Eq. 1), were used instead.We averaged 15-20 arterial voxels from the middle (MCA), posterior (PCA), and anterior cerebral arteries (ACA).We hypothesize that the resulting signal time course reflects vasodilation (given the negative signal change) at the arterial level (given the short delay and anatomical localization) in response to the CO 2 stimulus (see Results and Discussion for a more elaborate discussion of the AIF).
In DSC MRI, it is also useful to define an output function from the cerebral draining veins, whose magnitude can be used to scale the measured AIF 61 .Using both GRE-EPI and GRE-EPI-registered anatomical data, 15-20 venous voxels within or adjacent to the superior sagittal sinus (SSS), with S values above the 99.9th percentile, long delay (~ 2-3 s), and low S 0 (as vein contains more dOHb than artery, reducing the baseline GRE signal) were defined and averaged.This averaged venous time course is hereafter termed the venous output function (VOF).

Perfusion measurement
All voxel time courses were truncated to commence at the start of the AIF bolus and finish at the end of the VOF bolus.The AIF was then scaled by the integral of the VOF using the following equation: This step was conducted as, like in tissue, the VOF integral is representative of the amount of contrast agent (dOHb), whereas unlike in tissue, the AIF integral is representative of the degree of vasodilation.Normalizing to venous signal also reduces the effect of various perfusion quantification dependencies, such as the baseline oxygenation and susceptibility change induced by the contrast agent 11 .
Cerebral blood volume (CBV), cerebral blood flow (CBF), and mean transit time (MTT) maps were then calculated using a standard, truncated singular value decomposition (SVD) analysis, with an SVD noise threshold of 20%, brain density factor of 1.05, and hematocrit correction factor of 1.45 11,[30][31][32][33] .Although we employ an analysis based on a standard tracer kinetic model, it might be ideal for future work to incorporate ΔCBF into the analysis, particularly for the accurate estimation of baseline CBF and MTT.Note that the perfusion values are relative, not absolute; research has shown perfusion imaging values using DSC to depend on various contrast agent and acquisition parameters 10,11 .Thus, we report values in arbitrary units (a.u.).Also, note that all voxels yielding a negative CBV or CBF were removed from subsequent analysis, and are not included in the summary results (Table 1).

Delay measurement
Average bolus delay maps were also generated from the bhDSC data at 3 T and 7 T. Unlike the other perfusion metrics where DSC measurements were conducted in native space, the delay maps were calculated using the group average data.
First, the relaxation time courses were transformed to MNI152 2 mm anatomical space (using the previously described warp matrices) and averaged across subjects.The average time courses were then linearly interpolated to a temporal resolution of 0.5 s.The AIF time courses from each subject were then averaged.The subjectaveraged AIF was then shifted forward in 0.5 s steps, up to a maximum of 8 s, and the correlation between the AIF and tissue time courses for each voxel at each temporal shift was recorded.The temporal shift resulting in the highest correlation was determined to be the delay in a model-independent manner.A similar delay method has been employed in previous studies using resting-state fMRI data 78,79 .
Note that the calculated delay maps reflect a combination of artery-to-tissue delay, tissue transit time, and dispersion, given that the AIF time courses were simply shifted to find the maximal correlation with each voxel

Figure 1 .
Figure 1.Summary of the Breath-hold Experiment.(A) Description of the cerebrovascular reactivity phenomenon.At baseline, tissue blood oxygenation is between 60 and 80%.Increased CO 2 from breathholding causes vasodilation in the upstream arterial/arteriolar vasculature, leading to increased blood flow, and subsequently, a decrease in the downstream tissue dOHb concentration (hyperoxia), resulting in a GRE-MRI signal increase.(B) Signal time course (blue) resulting from nine consecutive breath-holds (green) in a representative venous voxel.(C) The signal time course resulting from the temporal averaging of the first eight boluses was then used for perfusion estimation.

Figure 2 .
Figure 2. Breath-Hold-Induced Relaxation Rate Time Course Dynamics.Time courses were averaged across subjects for gray matter (GM), white matter (WM), vein (VOF), and artery (AIF) at 3 T and 7 T.

Figure 3 .
Figure 3. Breath-hold CNR Properties.(A) CNR maps were calculated (Eq.2), transformed to MNI152 2 mm anatomical space, and averaged across subjects for both 3 T and 7 T (axial view is slightly dorsal to the lateral ventricles).For illustration purposes, CNR maps are shown when there is no bolus averaging (1) or 2, 4, 6, or 8 boluses averaged at 7 T. (B) Plot of subject-averaged CNR GM at 3 T (dashed) and 7 T (solid) as a function of the number of breath-holds averaged (error bars represent the standard error).Data are fit to a radical function ( a • √ b • x + c ) to demonstrate agreement with the well-known relationship between √ Averages and CNR.(C) Linear regression of subject-wise GM (red) and WM (blue) CNR values at 3 T vs 7 T. m represents the regression slope.

Figure 4 .
Figure 4.The Novel AIF.(A) Apparent ∆R 2 * AUC maps were calculated subject-wise, transformed to MNI152 2 mm anatomical space, and averaged at both 3 T and 7 T (sagittal and axial views are shown at the level of the middle cerebral artery).Pink arrows indicate an example of the AIF's location.AUC maps are overlaid onto the MNI152 2 mm anatomical template.(B) The associated apparent AIF ∆R 2 * time courses are shown at 3 T and 7 T, with subject-wise time courses in gray and subject-averaged time courses in pink.(C) Illustration of vasodilation mechanism, which results in arterial signal decrease (relaxation increase) during breath-holding.

Figure 5 .
Figure 5. bhDSC Perfusion Results.(A) CBV and CBF maps were calculated subject-wise, transformed to MNI152 2 mm anatomical space, and averaged at both 3 T and 7 T (sagittal (x = 30) and axial views (z = 38, 46) are shown).(B) Box and whisker plots show the calculated CBV and CBF values (for GM and WM) at 3 T and 7 T. (C) Voxel-wise regressions for CBV and CBF at 3 T vs 7 T.The regressions are based on voxels from the subject-averaged perfusion maps, and only voxels from the GM and WM are displayed.Regressions are set to intersect with the origin.m represents the regression slope.

Figure 6 .
Figure 6.bhDSC vs ASL.(A) CNR maps were calculated (Eq.2), transformed to MNI152 2 mm anatomical space, and averaged across subjects.(B) Box and whisker plots show calculated CNR values in GM for ASL and bhDSC at 3 T and 7 T. (C) CBF maps were calculated subject-wise, transformed to MNI152 2 mm anatomical space, and averaged across subjects.Axial view is slightly dorsal to the lateral ventricles for CNR and CBF maps.(D)Voxel-wise CBF regression between ASL and 7 T bhDSC.Regression is based on GM and WM voxels from the subject-averaged perfusion maps.As ASL coverage was variable between subjects, masking was conducted to only include voxels where ASL coverage overlapped amongst all subjects-refer to FigureS6for mask.We constrained analysis to axial slices 48 through 58, where ASL image intensity was homogeneous across each axial slice.Regression is set to intersect with the origin.

Table 1 .
Summary of bhDSC Perfusion Statistics.Average (mean), standard deviation (std), and coefficient of variation (CV) values for 3 T and 7 T are shown.