Cortical oxygen extraction fraction using quantitative BOLD MRI and cerebral blood flow during vasodilation

Introduction: We aimed to demonstrate non-invasive measurements of regional oxygen extraction fraction (OEF) from quantitative BOLD MRI modeling at baseline and after pharmacological vasodilation. We hypothesized that OEF decreases in response to vasodilation with acetazolamide (ACZ) in healthy conditions, reflecting compensation in regions with increased cerebral blood flow (CBF), while cerebral metabolic rate of oxygen (CMRO2) remained unchanged. We also aimed to assess the relationship between OEF and perfusion in the default mode network (DMN) regions that have shown associations with vascular risk factors and cerebrovascular reactivity in different neurological conditions. Material and methods: Eight healthy subjects (47 ± 13 years, 6 female) were scanned on a 3 T scanner with a 32-channel head coil before and after administration of 15 mg/kg ACZ as a pharmacological vasodilator. The MR imaging acquisition protocols included: 1) A Gradient Echo Slice Excitation Profile Imaging Asymmetric Spin Echo scan to quantify OEF, deoxygenated blood volume, and reversible transverse relaxation rate (R2 ’) and 2) a multi-post labeling delay arterial spin labeling scan to measure CBF. To assess changes in each parameter due to vasodilation, two-way t-tests were performed for all pairs (baseline versus vasodilation) in the DMN brain regions with Bonferroni correction for multiple comparisons. The relationships between CBF versus OEF and CBF versus R2’ were analyzed and compared across DMN regions using linear, mixed-effect models. Results: During vasodilation, CBF significantly increased in the medial frontal cortex ( P=0.004 ), posterior cingulate gyrus (pCG) ( P=0.004 ), precuneus cortex (PCun) ( P=0.004 ), and occipital pole ( P=0.001 ). Concurrently, a significant decrease in OEF was observed only in the pCG (8.8%, P=0.003 ) and PCun ( 8.7%,P=0.001 ). CMRO2 showed a trend of increased values after vasodilation, but these differences were not significant after correction for multiple comparisons . Although R2’ showed a slightly decreasing trend, no statistically significant changes were found in any regions in response to ACZ. The CBF response to ACZ exhibited a stronger negative correlation with OEF ( β=−0.104±0.027 ; t=−3.852,P<0.001 ), than with R2’ ( β=−0.016±0.006 ; t=−2.692,P=0.008 ). Conclusion: Quantitative BOLD modeling can reliably measure OEF across multiple physiological conditions and captures vascular changes with higher sensitivity than R2’ values. The inverse correlation between OEF and CBF across regions in DMN, suggests that these two measurements, in response to ACZ vasodilation, are reliable indicators of tissue health in this healthy cohort.


Introduction
The brain has a high metabolic demand for oxygen compared to other organs since the human brain comprises just 2% of the total body mass but consumes 20% of available oxygen for normal function (Rolfe and Brown, 1997), so regulation of cerebral blood flow (CBF) and oxygen delivery is critically important.Oxygen extraction fraction (OEF) is a key hemodynamic parameter to measure the brain's energy metabolism and is altered in aging and disease processes, including cerebrovascular disorders (Watchmaker et al., 2018;Fan et al., 2022;Liu et al., 2023), neurodegeneration (Jiang et al., 2020;Robb et al., 2022), and neuroinflammation (Cleland et al., 2021).In some pathological mechanisms, OEF is reduced concurrently with hypoperfusion due to decreased neural activity and lower oxygen metabolic demand (Ishii et al., 1996;Brown and Thore, 2011); while in other distinct mechanisms, OEF is elevated in areas of hypoperfusion, indicating physiological compensation in tissue that is at ischemic risk (Gupta et al., 2014;Fan et al., 2020;Lin et al., 2022).Therefore, it is important to measure OEF in tandem with CBF to understand overall oxygen consumption and brain tissue health.
In addition to resting state, imaging brain hemodynamics in varying physiological conditions, such as those created by controlled vasoactive stimuli, adds valuable information about cerebrovascular reactivity (CVR) and brain vascular health.These cerebrovascular "stress tests" monitor the CVR, or the brain's hemodynamic response to hypercapnic gas paradigms or to pharmacological vasodilation, e.g., with acetazolamide (ACZ) (Handwerker et al., 2007).CVR measurements reveal unique pathophysiological changes in intracranial stenosis (Lattanzi et al., 2018) and in aging that correlate with cognitive status (Silvestrini et al., 2006;Kim et al., 2021).During these vasoactive challenges, brain hemodynamics are typically monitored with the bloodoxygenation-level-dependent (BOLD) MRI signal due to its fast acquisition and ease of use; however, BOLD signal changes reflect multiple neuronal and vascular contributions and cannot be interpreted in the context of a specific hemodynamic parameter.Other investigations have used arterial spin labeling (ASL) MRI to directly assess the perfusion reactivity deficits in steno-occlusive disease and validate the BOLD response in the same patients (Mandell et al., 2008;Smeeing et al., 2016), highlighting ASL as a robust biomarker to quantify CBF in different physiological states and disease.Acetazolamide has been commonly used to assess the vasodilatory capacity of circulation because ACZ is more easily administered through intravenous or oral means than carbon dioxide, which requires a gas delivery system.In patients with major cerebral arterial occlusive diseases, ACZ also identified reduced reactivity in the hemisphere with the occlusive lesion even in patients who appeared to have preserved CO 2 reactivity (Kazumata et al., 1996).However, few studies have directly assessed OEF changes during vasoactive stimuli, largely because of a lack of non-invasive imaging approaches.[ 15 O]-PET scans have measured CBF changes in response to ACZ vasodilation in patients with cerebrovascular stenoses and showed a non-linear association between this perfusion reactivity and baseline OEF PET for various ischemic disease stages (Imaizumi et al., 2004;Nemoto et al., 2004;Yamauchi et al., 2004;Hokari et al., 2008).Unfortunately, these OEF PET measures are not easily available in different physiological states due to the complexity of the experiments and the need to administer multiple short-lived radiotracers for each condition.MRI-based OEF measures have recently been tested with ACZ administration and showed expected global OEF reductions after vasodilation using multiple T 2 relaxation-based measures (Baas et al., 2022).Venous blood oxygenation (Y v ), which can be derived from venous blood T 2 relaxation values in the sagittal sinus, combined with global cerebral metabolic rate of oxygen, has shown poorer oxygen utilization in patients with sickle cell disease compared to controls after vasodilation (Václavů et al., 2020), albeit with some bias between various T 2 -based MRI methods.Decreased OEF after ACZ vasodilation has also been demonstrated using 3-dimensional quantitative susceptibility mapping MRI (Buch et al., 2017).However, these values were averaged from larger, resolvable internal cerebral veins, thus limiting regional information on OEF response across the brain.
As an alternative to global OEF measurements using T 2 relaxation-based approach, quantitative blood-oxygenation-leveldependent (qBOLD) is a non-invasive MRI approach that models the temporal signal due to reversible transverse relaxation rate (R 2 ') to non-invasively extract OEF and deoxygenated blood volume (DBV) (Yablonskiy and Haacke, 1994).This method allows for voxelwise estimation of tissue OEF from asymmetric spin echo (ASE) scans, which can be useful in identifying areas of the brain that may be at risk for ischemia in stroke (Stone et al., 2019) or sickle cell disease (Guilliams et al., 2018;Wang et al., 2021).These clinical studies used qBOLD to identify regional, individualized pathophysiology in different ischemic tissue types (core, infarct growth, and contralateral tissue) in acute stroke.On the other hand, in patients with sickle cell anemia, abnormal OEF elevations were observed specifically in the deep white matter, consistent with microstructural damage, that were reduced postblood transfusion.Additionally, OEF activation has been reported to decrease in brain areas relating to motor task execution (Yin et al., 2018), especially under hypoxic conditions (Yin et al., 2021).In aging, studies have identified CVR reductions in brain regions that overlap with the default mode network (DMN), which is related to underlying vascular risk factors and suggests CVR in the DMN-CVR to be an important marker for brain health (Haight et al., 2015).Additional investigation of the reliability of qBOLD MRI in sensitive brain areas, including the DMN, across physiological states is thus critical to advance regional OEF as an informative biomarker in cerebrovascular disease and aging.
This study aimed to demonstrate a quantitative BOLD modeling MRI method for measuring OEF values in healthy participants during a physiological perturbation with ACZ.Previous studies have shown that OEF decreases in response to an increase in CBF from various vasodilating challenges, while CMRO 2 remains unaffected using different techniques (Kety and Schmidt, 1948;Vorstrup et al., 1984;Václavů et al., 2020;Vestergaard et al., 2023).We aimed to demonstrate this method by showing that OEF measured with qBOLD MRI is directly related to regional alterations in CBF.A secondary aim was to compare regional OEF with R 2 ' values in detecting an alteration with vasodilation and their correlation with perfusion.Furthermore, we focused on assessing the OEF relationship to perfusion in the DMN due to its functional associations with vascular risk factors (Tchistiakova et al., 2015) and cerebrovascular reactivity in cognitive impairment (Richiardi et al., 2015) to identify which DMN regions show the most significant OEF effect with vasodilation.

Study population
Eight healthy subjects (47 ± 13 years, 6 female) were recruited and gave written informed consent to participate in the study.
Participants must be at least 18 years old to participate in this study.Exclusion criteria included pregnancy, anemia, history of renal disease, hypertension, diabetes, stroke, or other known neurological diseases.Subjects with peripheral vasculopathy or Raynaud's disease that precludes IV administration and a history of allergy to sulfa drugs were excluded.In addition, participants that could not receive MRI due to the inability to lie motionless in the scanner, pacemakers, aneurysm clips, neurostimulators, artificial heart valves, metal objects in eyes, ear implants, and claustrophobia were excluded.This study was performed with approval from the Institutional Review Board at the University of California, Davis.

Vasodilatory stimulus
After baseline MRI scanning, all participants were given a slow administration dosage of 15 mg/kg ACZ as a pharmacological vasodilator over 2 min.After at least 15 min of uptake time to capture the full effect of pharmacological vasodilation, the participant received repeated scans in the post-ACZ state.
For OEF measurements, the Gradient Echo Slice Excitation Profile Imaging (GESEPI) ASE images were acquired with an echo planar imaging (EPI) readout with specialized z-shim gradients mitigating through-slice inhomogeneities (Blockley and Stone, 2016).Scan parameters included TR = 3,000 ms; TE = 56 ms; slice thickness = 1.25 mm; resolution = 2.3 mm 3 × 2.3 mm 3 × 1.25 mm 3 ; flip angle = 90 °; bandwidth = 2,004 Hz/Pixel; matrix size = 96 × 96 with 80 slices; field of view (FOV) = 224 mm with 7 values of the spin echo displacement time (τ) from 16 to 40 ms in the step of 4 ms.A total of 20 slabs were acquired, and each slab was constructed by averaging four 1.25-mm slices to correct for macroscopic field gradients with 100% partition oversampling (total 8 k-space partitions).The total acquisition time of the main ASE sequence was 2 min and 54 s to cover the whole brain with 80 slices.A separate spin-echo ASE scan (τ = 0) was collected with the same scan parameters as described, with an acquisition time of 24 s.
To generate quantitative CBF maps, we performed multi post-labeling-delay (PLD) 3D pseudo-continuous ASL (pcASL) MRI.The parameters were five different PLDs (0.9 s, 1.2 s, 1.4 s, 1.8 s, and 2.1 s) and effective labeling durations = 2 s; TR/TE = 4,700 ms/14 ms; slice thickness = 4 mm; spatial resolution = 3.4 mm 3 × 3.4 mm 3 × 4.0 mm 3 ; FOV = 220 mm (Li et al., 2015).Multi-band EPI readout was used with sliceacceleration factor of 6, FOV shift factor of 3, and EPI factor of 64.Including calibration scans with similar imaging parameters but no labeling to estimate longitudinal magnetization (M 0 ), the total acquisition time was approximately 6 min.The ASE and pcASL scans were acquired twice-once at baseline and once after vasodilation.

Image preprocessing
Prior to qBOLD modeling, ASE data were preprocessed first with motion correction of all echoes to the spin echo volume using MCFLIRT (Jenkinson et al., 2002).Next, the spin echo image was brain extracted using BET (Smith, 2002) to create a binary mask of the brain tissue, which was applied to all ASE R 2 '-weighted images.ASE data were finally smoothed with a Gaussian kernel (σ = 4 mm) in FSL (Jenkinson et al., 2012) to reduce the effect of noisy voxels.The empirical tradeoff between spatial resolution and noise reduction is shown in Supplementary Figure S4.
For pcASL data, images from all five PLDs were merged and underwent the same preprocessing steps as the ASE data.First, motion correction was applied with reference to the first pcASL image, followed by smoothing with a Gaussian kernel of 1.5 mm.Then, the M 0 image was brain extracted using BET to create the binary mask of the brain tissue, which was used for registration.

OEF quantification 2.4.2.1 Two-compartment model theory for quantitative BOLD
Quantitative measurements of the BOLD signal were used to non-invasively map hemodynamic parameters relating to brain metabolism and function.The quantitative BOLD (qBOLD) model (Yablonskiy and Haacke, 1994) measures oxygen extraction fraction from the measured reversible transverse relaxation rate, R 2 ' (where R 2 ' = R 2 *-R 2 ).A complete qBOLD model with R 2 '-weighted measurements can be achieved using asymmetric spin echoes (An and Lin, 2003), such as the GESEPI-ASE pulse sequence described above.qBOLD modeling was performed with two compartments (tissue and blood), as this approach provides reliable estimates of OEF with a full range of echoes for fitting, including those before the spin echo (Cherukara et al., 2019).
The objective of qBOLD modeling is to separate OEF and DBV effects from the ASE signal and generate a brain map of OEF.The (1) tissue signal, (Yablonskiy and Haacke, 1994), S t is: where δω is the characteristic frequency (with R 2 ′ DBV × δω), t c is the characteristic time, and R t 2 is the irreversible transverse relaxation rate of bulk tissue.The R 2 '-weighted signal has different behavior in two different regimes of the parameter τ (Eq.1).The boundary between these regimes is considered to be at T c = 1.76t c (Cherukara et al., 2019).The characteristic time t c was defined as t c 1 δω (Lee et al., 2018).The time variable τ corresponds to the spin echo displacement values of the ASE acquisition.The (2) blood signal, S b , (Sukstanskii and Yablonskiy, 2001;Yablonskiy et al., 2013) is: where C(η) and S(η) are Fresnel function, and η ( 3δω|τ| π ) 2 .R b 2 is the transverse relaxation rate of blood and is described as a function of the fractional hematocrit (Hct) and OEF of the intravascular compartment (Simon et al., 2016): The total signal calculated from a voxel in this two-compartment tissue model is the sum of the signal from each compartment: After fitting the two-compartment tissue qBOLD model, R 2 ' and DBV were estimated.Then, OEF was determined by the relationship between R 2 ' and DBV with known constants: where Hct is the patient's fractional hematocrit, Δχ 0 is the magnetic susceptibility difference between oxygenated and deoxygenated red blood cells, B 0 is the magnetic field, γ is proton gyromagnetic ratio.R 2 '-weighting was acquired by shifting the spin echo refocusing pulse for 7 values of the spin echo displacement time, τ, in the step of 4 ms from 16 ms to 40 ms.The spin echo (τ = 0 ms) was collected in separate scan.Hct was assumed at 0.40 (Nicoll et al., 2012).Other physiological parameters were also set as constant values as follows: Δχ 0 0.264 × 10 −6 (Spees et al., 2001); R t 2 11.5s −1 (He and Yablonskiy, 2007).R t 2 11.5s −1 .Additionally, acquisition parameter values used for the model were B 0 = 3T, TE 56ms, and TR 3s.

Bayesian inference for mapping OEF with qBOLD
Fitting the acquired ASE signals to the two-compartment tissue model allows us to separate OEF (our parameter of interest) and DBV per voxel.To fit our ASE signal data, we used a Bayesian framework that estimates R 2 ' and DBV, which have shown to be more separable parameters with a suitable posterior distribution for the estimation (Cherukara et al., 2019).Given the total signal S Total in Eq. 4, the Bayesian approach will find the optimal pair of (DBV and R 2 ') using a variational Bayesian (VB) inference scheme implemented in FAST ASL and BOLD Bayesian Estimation Routine (FABBER) toolkit (Woolrich et al., 2006;Chappell et al., 2009;Groves et al., 2009) in FSL (version 6.0.5, Oxford, United Kingdom), as illustrated in Figure 1A.VB inference was performed with fixed prior mean values, μ 0 , and standard deviations, σ 0 The mean values for R 2 ' and DBV for healthy subjects were 2.6 s −1 and 3.6%, respectively (Stone and Blockley, 2017).The optimal prior standard deviation was chosen as σ 0 (R 2 ′ ) 10 3 2 s −1 and σ 0 (DBV) 10 1 2 % based on numerical simulations (Le et al., 2023).After implementing Bayesian estimation to map DBV and R 2 ' in each participant, OEF was calculated voxelwise based on the proportion between DBV and R 2 ' with known constants from Eq. 5.All qBOLD model analyses were done in native space.

The relationship between OEF and CMRO 2
The relationship between CMRO 2 (μmol/100 g/min) and OEF (%) (Kety and Schmidt, 1948) can be expressed as Where CBF is the cerebral blood flow (mL/100 g/min), and [Hb] a is the oxygenated hemoglobin concentration in the arteriole (8.272 μmol/mL) This is calculated from [Hb] a [Hb] t × Y a , where [Hb] t (8.441 μmol/mL) is the hemoglobin concentration in the tissue blood with assumption of hematocrit of Hct = 0.40 (Zhang et al., 2018).Y a is the arterial oxygenation, which is assumed to be 0.98.

Perfusion measurements
The complete procedure to measure cerebral blood flow is shown in Figure 1B.We used the Bayesian Inference for Arterial Spin Labeling MRI (BASIL) toolkit (Chappell et al., 2009) implemented in FSL (version 6.0.5, Oxford, United Kingdom) for analyzing ASL data and quantifying perfusion maps.This approach allows fitting a two-compartment model to separate macrovascular and tissue information in ASL signal to quantify CBF maps (Groves et al., 2009;Chappell et al., 2010).
The two-compartment model assumed that the walls of arterial vessels do not allow any substance to pass through them, and that blood in the arteries moves through the voxel immediately.The total signal from any voxel is the sum of both tissue and intravascular components.CBF maps underwent a voxelwise calibration method using the measured M 0 maps (Pinto et al., 2020).Fix label duration was applied, and initial parameter values were assumed as following: initial arterial transit time (ATT) = 1.3 s; T 1 values for tissue (T 1t ) = 1.3 s (based on 3 T field strength); T 1 values for blood (T 1b ) = 1.65 s; brain/blood partition coefficient (λ) = 0.9 mL/g; and labeling efficiency (α) = 0.85 (Alsop et al., 2015).Analyses to compute CBF maps were done in native ASL space.

Image registration
To assess the regional relationship between hemodynamic parameters, especially OEF and perfusion across subjects, all quantitative maps were registered to Montreal Neurological Institute (MNI) space (Collins et al., 1999;Fonov et al., 2009;Fonov et al., 2011).For hemodynamics parameter maps, including OEF, DBV, and R 2 ', the spin echo image from the ASE scans was linearly co-registered, with FMRIB's Linear Image Registration Tool (FLIRT) (Jenkinson and Smith, 2001;Jenkinson et al., 2002;Greve and Fischl, 2009), to subjectspecific T 1 -structural space and then non-linearly coregistered, with FMRIB's Non-linear Image Registration Tool (FNIRT) (Andersson et al., 2010), to MNI space.All parameter maps were then warped to MNI space with combined transformation parameters from those two registrations.To register CBF maps to MNI space, we first performed linear registration (FLIRT) of proton density (M 0 ) images to T 1 images; then, all CBF maps were transformed into MNI space using non-linear co-registration as described above.

Mean ROI calculation
Before assessing parameters of interest (CBF, OEF, and R 2 ') across brain regions and calculating the average for each ROI, we excluded voxels with unphysiological extreme values.Voxels with OEF values greater than 100%, and with estimated R 2 ' greater than 20 s −1 were excluded from analysis.The thresholds were set consistently to prior relevant studies (Stone and Blockley, 2017;Cherukara et al., 2019).For group average parameter maps across all participants, we calculated the mean only for MNI voxels having at least 50% of the population (i.e., at least 4 subjects) within physiological values.Intermediate processing steps, including subject masks in MNI space, are shown in Supplementary Material, Supplementary Figure S1 (for OEF), Supplementary Figure S2 (for DBV), and Supplementary Figure S3 (for R 2 ').

Statistical analysis
All statistical analyses were performed using MATLAB R2021a (MathWorks, Natick, MA, 2016).Differences between before and after vasodilation were assessed using a pairwise two-way t-test with Bonferroni correction for nine ROIs.A p-value of less than 0.006 (where α 0.05 number of comparisons ) was considered statistically significant for this analysis.
To evaluate the relationship of OEF and CBF across regions during vasodilation, a linear mixed-effects model was conducted using OEF as the dependent variable and CBF as the independent variable, with subject and region as random effects.We also investigated the correlation of R 2 ' and CBF across regions during vasodilation, using a separate linear mixed-effects model with R 2 ' as the dependent variable and CBF as the independent variable.For both models, subjects and regions were used as random effects to adjust for inherent physiological variations across individuals and regions.Each mixed model was fit separately using the MATLAB software package "fitlme" (Pinheiro and Bates, 1996).A p-value of less than 0.05 was considered statistically significant.

Results
All eight healthy participants received MRI scans consisting of the GESEPI-ASE sequence for quantitative BOLD modeling to estimate hemodynamic parameters, including R 2 ', DBV, and OEF (our parameter of interest).After excluding unphysiological voxels, group average maps of each parameter before and after vasodilation are shown in Figure 2. Excluded voxels were primarily located in frontal brain areas prone to air-tissue bulk magnetic susceptibility effects.The percentage of voxels removed due to thresholding in the two conditions is shown in Supplementary Table S1 (Pre-ACZ: DBV = 0.21%, R 2 ' = 2.05%, OEF = 13.1%;Post-ACZ: DBV = 0.16%, R 2 ' = 1.84%,OEF = 10.93%.We focused on parameters that are sensitive to brain oxygenation (OEF and R 2 ') and perfusion; the corresponding DBV maps are shown in Figure 2.
Compared to the baseline condition, increased CBF and decreased OEF were observed across multiple areas in the brain at the group average level (Figures 2A, B).We also observed a decreasing trend in R 2 ' (Figure 2C) and an increasing trend in DBV (Figure 2D) during vasodilation.Mean parameter values, including CBF, OEF, R 2 ', DBV and CMRO 2 for all subjects across selected regions for each condition are presented in Table 1.Additionally, the distribution of changes in OEF, R 2 ', DBV and CBF in response to ACZ for each selected region were plotted in Figure 3.
Paired two-way t-test with Bonferroni correction showed statistically significant differences in CBF and OEF estimates in selected regions during vasodilation (i.e., Pre-ACZ versus Post-ACZ), as presented in Table 1.With the injection of ACZ, CBF increased across all chosen regions before correction for multiple comparisons.After correction, elevated CBF remained significant in MFG (P 0.004), pCG (P 0.004), PCun (P 0.004), and OP (P 0.001) regions.In parallel, a significant reduction in OEF from qBOLD was observed only in pCG (P 0.003) and PCun (P 0.001) after correction, with a trend of decreased OEF also in the angular gyrus.In addition, a slight trend of DBV increase during vasodilation was observed only in pCG (P 0.045) before correction.For R 2 ', there was a small reduction trend during vasodilation; however, we found no region showing statistically significant change with ACZ (Table 1).Additionally, while there was a trend of increased CMRO 2 after vasodilation in several brain areas, especially the MFG and MTG, CMRO 2 showed no significant change in all regions after correction for multiple comparisons (Table 1).
The CBF response to ACZ showed a significant negative linear relationship to OEF value at baseline and after vasodilation using a mixed-linear effects model with subject clustering (β −0.104 ± 0.027; t −3.852, P < 0.001) (Figure 4A; Table 2).We further investigated the relationship between CBF and R 2 ' using the same mixed-effects model.The CBF response to ACZ was also correlated with R 2 ' (β −0.016 ± 0.006; t −2.692, P 0.008) (Figure 4B; Table 2).

Discussion
The purpose of this study was to demonstrate qBOLD modeling of ASE scans in healthy volunteers to detect regional OEF decreases in key ROIs of the DMN during concomitant perfusion increase with ACZ vasodilation.The main findings were as follows: 1) Across baseline and vasodilation states, OEF was inversely related to quantitative perfusion as expected, indicating a compensation to maintain oxygen metabolism and reliability of the local OEF measures.2) Similar changes with ACZ observed using R 2 ' relaxation estimates from the same ASE acquisitions, which highlights the utility of qBOLD compartment modeling to non-invasively measure specific OEF imaging markers.3) OEF reduction was observed in critical DMN regions, including the precuneus and posterior cingulate gyrus.

Baseline OEF and change with acetazolamide in healthy participants
The observed mean baseline OEF across selected brain regions of 35.1% is consistent with the physiological range of previous PET and MRI studies (Yamaguchi et al., 1986;Leenders et al., 1990;Cho et al., 2021;Jiang et al., 2023).We observed an OEF reduction of 6.2%-8.7%(absolute oxygenation) in significant DMN regions, which was smaller in magnitude than the 15.6% absolute OEF decrease in the sagittal sinus from past ACZ studies on healthy controls (Václavů et al., 2020).This discrepancy could be due to methodological differences, as previous studies focused on global OEF values and used a distinct contrast mechanism and calibration with T 2prepared inversion recovery sequences (Václavů et al., 2020;Baas et al., 2022).The qBOLD method adopted here may also have lower sensitivity to OEF changes due to noise contributions in the model fitting in vivo.Of note, the T 2 -based study also reported a larger magnitude of perfusion increase (69.3% increase) using pcASL with ACZ than our study (average 43.4% increase across ) (in s −1 ), deoxygenated blood volume (DBV) (in %) and cerebral metabolic rate of oxygen (CMRO 2 ) (in μmol/100 g/min) before (Pre) and after (Post) vasodilation across all healthy subjects (mean ± std) in different regions of interest (ROIs): angular gyrus (AG), medial frontal gyrus (MFG), anterior cingulate gyrus (aCG), posterior cingulate gyrus (pCG), precuneus (PCun), occipital pole (OP), supramarginal gyrus (SG), middle temporal gyrus (MTG), and inferior temporal gyrus (ITG).All tests were performed using a two-sided paired t-test with Bonferroni correction.(*p < 0.05, significant raw p-value; **p < 0.006, significant p-value after correction).regions).Therefore, direct comparison of the OEF imaging approaches in the same participants and improvement to sensitivity of our qBOLD approach is warranted in future work.

Correlation between hemodynamic parameters and perfusion
Although R 2 ' has been used as a surrogate marker for brain OEF in previous studies (Ni et al., 2015), including identification of putative oxygenation changes in cerebrovascular disease (Ni et al., 2017;Kaczmarz et al., 2021), this work suggests that qBOLD compartmental modeling outperforms R 2 ' measures in detecting physiological OEF changes in healthy volunteers.While significant OEF reductions were measured in the pCG and PCun with qBOLD, none of the ROIs showed a significant R 2 ' change with ACZ.Additionally, the expected inverse relationship between OEF and CBF across baseline and vasodilatory states was stronger using qBOLD than between R 2 ' and CBF.One rationale could be that R 2 ' is proportional to both OEF and DBV, such that OEF reductions concurrently with potential DBV elevations during vasodilation lead to overall smaller changes in R 2 ' values.Because qBOLD compartment modeling uses multiple asymmetric echoes to disentangle contributions from OEF and DBV, the resulting regional OEF measures are likely to be more sensitive than R 2 ', as observed in our results.In addition to the Bayesian fitting for R 2 ', a supplementary analysis of the monoexponential fits based on streamline qBOLD (Stone and Blockley, 2017) was also performed.Our more complex modeling approach may produce R 2 ' measures that are noisier compared to monoexponential fits for the relaxation parameter (Supplementary Figure S5), which may contribute to lack of sensitivity of R 2 ' to changes during vasodilation.Besides, in whole group level analysis with monoexponential fits, there was a slight decrease for R 2 ' in pCG (P 0.028) and OP (P 0.039) before correction (Supplementary Table S2).However, in this analysis, we chose to compare OEF and R 2 ' values that were both derived from Bayesian fitting for consistency.R 2 ' measures are also susceptible to contributions from orientation effects of the blood vessels (Kaczmarz et al., 2021) and artifacts from non-heme susceptibility sources (iron, myelin).These contributions will also manifest differently on OEF values from qBOLD modeling, which requires further evaluation with related multi-parametric qBOLD approaches that acquire separate R 2 and R 2 * maps measures to assess oxygenation (Gersing et al., 2015).Such multi-parametric qBOLD have been used to identify reduced baseline oxygen metabolism in the affected middle cerebral artery territory with moderate ischemia and may provide an alternative to ASE acquisitions (Bouvier et al., 2015).
A good validation for our qBOLD method to measure OEF is to test the hypothesis that CMRO 2 is unchanged during vasodilation, i.e., that OEF decreases are commensurate with the increase in CBF during acetazolamide.From this study, we found no significant CMRO 2 change after Bonferroni correction in all selected regions, which is consistent with our hypothesis and reflects reliability of the qBOLD method.However, we observed a slight increasing trend in CMRO 2 , which may reflect residual underestimation bias in OEF values (i.e., underestimating the OEF decrease) and limited sensitivity of the qBOLD model to OEF.While CBF significantly increased in all selected regions, OEF only showed a significant decrease in pCG and PCun (Table 1).Therefore, there was an overall trend of increased CMRO 2 in several regions, even though we observed the expected compensation of OEF and CBF during vasodilation.Future work will enhance the accuracy and sensitivity of qBOLD modeling with refined physiological prior information and improved fitting routines that are more robust to noise contributions.

Effects of vasodilation in cerebral regions
This study imaged the regional effects of vasodilation on OEF across the entire cerebral cortex, focusing on DMN (Figures 2, 3).We found that most DMN regions showed increased CBF during vasodilation, and decreased OEF was observed in pCG and PCun.In contrast, the non-DMN regions did not show a significant response in any hemodynamic parameters.One potential implication of this finding is that DMN regions show a strong sensitivity to vasodilation in those regions, which is consistent with well-established mechanisms of vasodilation with ACZ (Kleinschmidt et al., 1995;Wang et al., 2015).In addition, BOLD fMRI studies during hypercapnia have identified co-fluctuations in brain hemodynamics that spatially resemble the DMN but are responsive to vasoactive (not neuronal) stimuli (Bright et al., 2020).Such "vascular networks" may have unique OEF properties at baseline and in vasodilated conditions, and methods such as the qBOLD approach can characterize such physiological mechanisms of brain functional network emergence.These capabilities will enable important future studies to assess the link between oxygen metabolism in other brain networks and their disorders.Outside of the DMN, the medial temporal lobe (MTL) has been known as a key early region that is affected by many brain diseases (Clifford R. Jack et al., 1997;Tam et al., 2005;Burton et al., 2009).OEF can serve as a functional biomarker due to its sensitivity in the early stages of brain diseases and aging, as neural activity is tightly linked to the brain's oxygen consumption (Watts et al., 2018).At baseline, average OEF values in the medial temporal lobe were measured to be 25.9% ± 3.0% (Table 1), consistent with prior MRI studies' measurement of OEF in vessels supplying the MTL of 23.9% using T 2 -relaxation-under-phase-contrast approach (Jiang et al., 2023).Even though there was a decreasing trend, we observed no significant change in MTL OEF during vasodilation (P 0.609).This observation contrasts with Jiang et al., who observed increased MTL OEF after vasoactive challenge with caffeine ingestion.Caffeine is a common vasoconstrictor, which is known to reduce CBF and increase OEF to maintain the same amount of total oxygen consumption (Xu et al., 2015).Both ACZ and caffeine are vasoactive stimuli with minimal neuroactive effects, although they are in opposite directions.For the MTL, Jiang et al. observed a 9.1% OEF increase due to caffeine, while our study observed an average of 12.7% OEF decrease across regions due to ACZ vasodilation.More studies are needed to characterize the accuracy of qBOLD to detect changes in MTL OEF with physiological changes or aging and may require additional corrections for signal loss or distortion due to susceptibility effects in some MTL subregions (Olman et al., 2009).

Limitations
To perform qBOLD modeling, we have used an ASE sequence to directly estimate R 2 ' with consistent R 2 -weighting across echoes (An and Lin, 2003).One main limitation of the qBOLD model is that the separation of DBV and OEF effects in the qBOLD model depends mainly on the subtle change in decay patterns.Therefore, it requires high signal-to-noise-ratio (SNR) to accurately estimate OEF from the ASE signals (An and Lin, 2000).Due to challenging low SNR from our sequence, OEF values estimated from qBOLD model still had unreasonable physiological values in particular voxels that we had to remove.To address this noise limitation, a benefit to fitting this multi-compartment with the Bayesian framework is the ability  to include prior knowledge to improve parameter estimates.In this study, prior mean values and standard deviations were taken from previous studies (Stone and Blockley, 2017;Cherukara et al., 2019;Le et al., 2023).Moreover, it has been suggested that prior means do not significantly affect the parameter estimation (Cherukara et al., 2019).However, the prior standard deviations were optimized based on a broader range of R 2 '-weighted images (24 R 2 '-images in Cherukara et al., 2019, and14 in Le et al., 2023) than acquired in this MRI protocol.For this data, we also investigated the effect of prior standard deviations to the parameter estimation by comparing between the chosen prior standard deviations (σ(DBV) 10 3 2 %; σ(R 2 ′ ) 10 1 2 s −1 ) with a broader one (σ(DBV) 10 5 2 %; σ(R 2 ′ ) 10 5 2 s −1 ).We found that the broader standard deviation, which places less weight on prior information, resulted in more unphysiological values compared to our chosen value as shown in Supplementary Table S1.Further optimization for prior information will be assessed to maintain a robust estimation of OEF.
Hct was assumed as 0.40 based on the literature value for general circulation (Nicoll et al., 2012), which is consistent with average hematocrit in healthy subjects from other studies (Brown et al., 1962;Jeong et al., 2021).On the other hand, prior relevant studies have used a lower Hct of 0.34 (He and Yablonskiy, 2007) because of lower Hct in small vessels.Because Hct influences quantification specifically of the intravascular compartment, for our twocompartment model, changing Hct has a non-linear effect on the parameter estimates.For instance, using a Hct of 0.34 (i.e., a 15% lower value of Hct) and without spatial regularization, Cherakura et al. observed the same mean DBV values but a 3% higher mean OEF values compared to analysis with Hct of 0.40.However, it is unclear how much the Hct varies throughout the brain, which may propagate to additional regional differences of OEF estimates from the qBOLD signal.Finally, this study adopted a Bayesian framework that aims for the maximal free energy of the approximate posterior to better estimate the data (Chappell et al., 2009).This approach was chosen to ensure that applying additional smoothing (Gaussian smoothing with kernel of 4 mm) did not affect the fitting.Because the qBOLD-OEF measurement requires more complex models, it is more sensitive to the presence of noise compared to CBF quantification from ASL.Therefore, we applied a larger kernel of smoothing (i.e., 4 mm) for OEF, which is consistent to previous studies (Cho et al., 2021;Uchida et al., 2022).This kernel was empirically chosen after comparing multiple kernel sizes on our data to reduce unstructured noise and heterogeneity in the OEF maps but has the drawback of decreasing the spatial resolution of regional OEF measures.The relatively small sample size of our study is also a limitation and should be expanded in future studies to include larger cohorts with various ages and cerebrovascular conditions.

Conclusion
We demonstrated that regional OEF measurements from qBOLD are reliable across multiple physiological conditions and evaluated OEF through correlation with ASL perfusion in response to ACZ vasodilation as reliable indicators of tissue health in healthy cohort.This study enhances our understanding of the relationship between OEF and CBF with regional alterations in brain vasculature and has implications for the future use of oxygenation MRI to study vascular alterations in cognitive impairment and cerebrovascular diseases.
FIGURE 4 (A) The correlation between oxygen extraction fraction OEF (%) and cerebral blood flow (CBF); (B) the relationship between R 2 ' (s −1 ) and CBF across all healthy subjects in nine ROIs (red: before vasodilation (Pre-ACZ); blue: after vasodilation (Post-ACZ).The solid black line represents the fitted line from mixed-effect models.Linear mixed-effect models with adjustment for subject clustering and region were performed separately for OEF and R 2 ' .p < 0.05* Mixed-effects model: OEF~CBF + (1|Region) + (1|Subject).

TABLE 2
Linear mixed-effects models between OEF, R 2 ', and CBF in the whole group (N = 8) across 9 ROIs.In both models, subjects and regions were used as random effects.*p < 0.05 β: standardized beta coefficient; SE: standardized error; OEF: oxygen extraction fraction; R 2 ': transverse relaxation rate; CBF: cerebral blood flow.