Lag-Optimized Blood Oxygenation Level Dependent Cerebrovascular Reactivity Estimates Derived From Breathing Task Data Have a Stronger Relationship With Baseline Cerebral Blood Flow

Cerebrovascular reactivity (CVR), an important indicator of cerebrovascular health, is commonly studied with the Blood Oxygenation Level Dependent functional MRI (BOLD-fMRI) response to a vasoactive stimulus. Theoretical and empirical evidence suggests that baseline cerebral blood flow (CBF) modulates BOLD signal amplitude and may influence BOLD-CVR estimates. We address how acquisition and modeling choices affect the relationship between baseline cerebral blood flow (bCBF) and BOLD-CVR: whether BOLD-CVR is modeled with the inclusion of a breathing task, and whether BOLD-CVR amplitudes are optimized for hemodynamic lag effects. We assessed between-subject correlations of average GM values and within-subject spatial correlations across cortical regions. Our results suggest that a breathing task addition to a resting-state acquisition, alongside lag-optimization within BOLD-CVR modeling, can improve BOLD-CVR correlations with bCBF, both between- and within-subjects, likely because these CVR estimates are more physiologically accurate. We report positive correlations between bCBF and BOLD-CVR, both between- and within-subjects. The physiological explanation of this positive correlation is unclear; research with larger samples and tightly controlled vasoactive stimuli is needed. Insights into what drives variability in BOLD-CVR measurements and related measurements of cerebrovascular function are particularly relevant when interpreting results in populations with altered vascular and/or metabolic baselines or impaired cerebrovascular reserve.


INTRODUCTION
Cerebrovascular Reactivity (CVR), the cerebral blood flow (CBF) response to vasoactive stimuli, reflects the regulatory ability and health of the cerebrovasculature (Battisti-Charbonney et al., 2011;Meng and Gelb, 2015). Using Blood Oxygenation Level Dependent functional MRI (BOLD-fMRI) to reflect CBF changes is common for CVR mapping (Pinto et al., 2021). However, theoretical and empirical evidence shows that task-induced and resting-state BOLD amplitude changes are modulated by the vascular and metabolic baseline state (Lu et al., 2008;Griffeth et al., 2011;Kim and Ogawa, 2012;Liu, 2013;Chu et al., 2018). Multiple studies have artificially decreased or increased a subject's resting CBF and then compared BOLD task activation amplitudes across these different states: a higher resting CBF corresponds to a lower task-based BOLD amplitude change (Cohen et al., 2002;Brown et al., 2003;Stefanovic et al., 2006;Vazquez et al., 2006). These studies primarily focus on BOLD signal changes evoked by neural activity; for BOLD-CVR effects evoked by altered blood gas levels, the literature gives less consensus on the directionality of this baseline modulation. However, there is still evidence that the baseline vascular state is coupled with CVR. For instance, positive correlations between baseline cerebral blood flow (bCBF) and BOLD-CVR have been reported across individuals (Leoni et al., 2017), and studies that include both bCBF and CVR show they both decrease with age throughout certain years of adulthood (Lu et al., 2011;Leung et al., 2016). In healthy adults, there is evidence that CVR amplitudes and timings can be modulated experimentally by changing the level of baseline vasodilation or vasoconstriction (Bright et al., 2011;Halani et al., 2015).
In pathological conditions, the capacity for a vessel to dilate in response to a stimulus can be diminished due to a predilated baseline, for example, in sickle cell disease bCBF is shown to be elevated and CVR diminished (Kosinski et al., 2017;Václavů et al., 2019). In steno-occlusive diseases, a vasodilatory stimulus may cause paradoxical decreases in blood flow to an area with exhausted dilatory reserve, as vascular resistance is reduced in surrounding regions (Sobczyk et al., 2014). This observation is known as the "vascular steal" phenomenon and can manifest as apparent negative CVR responses, with and without bCBF alterations (Sobczyk et al., 2014). The mechanistic relationship between bCBF and CVR in pathological cases is often not straightforward; to bring clearer interpretations to these pathological cases, a better understanding of the relationship between bCBF and BOLD-CVR in healthy populations is needed.
Variability in CVR in healthy young populations, modulated by genetic risk factors, could be a significant predicter of neurological pathology in later life (Suri et al., 2015). Therefore, understanding how baseline vascular physiology relates to dynamic cerebrovascular processes will bring further insights into what drives participant variability in BOLD-CVR measurements and related measurements of cerebrovascular function. This improved understanding will help to measure and track cerebrovascular health across the lifespan. Furthermore, BOLD-fMRI signal changes are commonly used to infer changes in neural activity, yet resting-state and task-based BOLD-fMRI amplitudes can be strongly influenced by CVR (Chen and Gauthier, 2021;Cohen et al., 2004;Liu et al., 2013a). Therefore, validity of using BOLD-fMRI as a surrogate measure of neural activity depends on the relative vascular and neural contributions to the BOLD-fMRI signal, and our ability to separate these contributions (Kannurpatti et al., 2010;Tsvetanov et al., 2020;Bright et al., 2020). These factors are particularly important to consider in clinical cohorts that present with altered vascular and/or metabolic baselines.
Here, we add to the small body of literature assessing bCBF and BOLD-CVR relationships, and address how key CVR data acquisition and modeling choices affect this relationship. Specifically, unlike most previous literature, we investigate the relationship between naturally varying bCBF and BOLD-CVR, in a healthy young population, and not in situations where bCBF has been artificially altered or in the presence of aging or pathology. We characterize both the between-subject relationship between bCBF and CVR in gray matter (GM), and the within-subject spatial relationship of bCBF and BOLD-CVR across different cortical brain areas.
Our previous work assessed a practical modification to a typical resting-state BOLD-fMRI protocol, showing that the addition of a short breathing task to a resting-state fMRI acquisition facilitated CVR modeling and quantification. The inclusion of a non-invasive breathing task (breath-hold or cued deep breathing (Pinto et al., 2021)) increased the number of voxels in the brain that had a significant relationship between the partial pressure of end tidal CO 2 (P ET CO 2 ) and BOLD-fMRI, and improved our confidence in these voxel-wise CVR estimates. The breathing tasks evoke large, transient changes in the arterial blood CO 2 pressure, which consequently changes blood flow; the P ET CO 2 measurement is our surrogate measure of arterial blood CO 2 changes (Takano et al., 2003;McSwain et al., 2010). If BOLD-CVR maps reflect sensible physiological variation we would hypothesize some spatially sensitive relationship with bCBF, based on previous work cited and discussed in the first paragraph of this introduction. In this paper, we assess the impact of methodological choices on between-and within-subject correlations of bCBF and BOLD-CVR. Specifically, we predict that adding breathing tasks to drive robust P ET CO 2 changes will improve our detection of relationships between BOLD-CVR and bCBF. Furthermore, to achieve accurate estimates of CVR amplitude, it is important to consider spatially variable hemodynamic delays (lags) between the P ET CO 2 regressor and the BOLD signal, by shifting the P ET CO 2 regressor in time to optimize the full model fit (Moia et al., 2020;Stickland et al., 2021). There is a temporal offset between the P ET CO 2 recording and the fMRI recording, and this is driven by both methodological factors (delay between CO 2 exhalation inside the scanner and the recording of P ET CO 2 in the control room) and physiological factors (brain-to-lung, vascular transit delays, variability in the vasodilatory capacity of local arterioles and the spatiotemporal complexities of the BOLD response). It is therefore important to estimate the CVR amplitude alongside a correction for variable CVR timing, on a regional basis (Stickland et al., 2021). Based on this, we predict that lag-optimization of BOLD-CVR estimates will lead to a stronger relationship between bCBF and BOLD-CVR. Note, this manuscript focuses on the relationship between CVR maps presented previously (Stickland et al., 2021) and newly presented bCBF maps.
A pseudo-continuous arterial spin labelling (pCASL) dataset was also acquired whilst the subject was at rest, with sequence parameters guided by the white paper (Alsop et al., 2015): a background-suppressed 3D GRASE read-out, 5 segments, TR/TE = 4000/19.4 ms, 4 mm isotropic voxels, 40 axial slices (interleaved, ascending), FOV = 256 mm, post label delay/Tau duration = 1800 ms, 90 mm label offset from middle slice. This sequence had 11 tag and 11 control volumes and one M0 volume, for a total scan time of approximately 7.5 min.
Inspired and expired CO 2 and O 2 pressures (in mmHg) were sampled through a nasal cannula worn by the participant, recorded at 1000 Hz with LabChart software (v8.1.13, ADInstruments), connected to a ML206 Gas Analyzer and PL3508 PowerLab 8/35 (ADInstruments) and synced with volume triggers from the scanner.

Data Analysis
The data from this study unfortunately cannot be made openly available due to restrictions of the ethical approval that they were collected under. However, analysis derivatives that are not included in this manuscript may be provided, on request, within ethical guidelines.

Blood Oxygenation Level Dependent Cerebrovascular Reactivity Analysis
For a detailed description, please refer to our previous work where we present the analysis of these CVR maps using the same subjects (Stickland et al., 2021). In brief, volume registration was performed on each fMRI dataset, with the SBRef volume from the middle (second) fMRI dataset used as the reference volume. Brain extraction was then performed on each fMRI dataset. P ET CO 2 values were identified, convolved with a canonical hemodynamic response function (hrf), and shifted ± 15 s in 0.3 s increments, then down-sampled to the TR. Multiple linear regression was performed separately for five data segments, two of which included breathing modulations and three of which only included resting-state data. These five data segments were created from the three acquisitions as illustrated in Figure 1, which also displays the P ET CO 2 trace and the GM average BOLD-fMRI trace for each acquisition. The regression model consisted of mean, drift terms, 6 motion parameters, and a P ET CO 2 time-series. The beta-weight for the unshifted P ET CO 2 regressor, scaled by the fitted mean, produced CVR maps with no lag optimization (No-Opt CVR, units:%BOLD/mmHg). For CVR maps with lag optimization (Lag-Opt CVR), the model was run for each shifted P ET CO 2 regressor; parameter estimates were taken from the model with the largest full model R-squared (Moia et al., 2020;Stickland et al., 2021).

Baseline Cerebral Blood Flow Analysis
First, volume registration was run with AFNI's (Cox, 1996) 3dvolreg, and then the M0 volume was separated from the tagcontrol volumes. Brain extraction with FSL was performed on the M0 image to create a brain mask. FSL's BASIL toolbox (Chappell et al., 2008) was used for perfusion modeling and quantification of voxel-wise cerebral blood flow in ml/100 g/min. Guided by the ASL white paper (Alsop et al., 2015) the analysis inputs included: inflow time of 3.6 seconds (post label delay of 1.8), bolus duration of 1.8 seconds, T1 of 1.3 seconds, adaptive spatial smoothing, inversion efficiency of 0.85, M0 brain masking and voxel-wise calibration with the M0 image. The fsl_anat directory was given as an input to BASIL, which performed registration to T1-weighted space.

Transformation of Atlas and Cerebrovascular Reactivity Parameter Files to T1-Weighted Space
In FSL, left and right hemisphere brain masks were created based on the MNI152 6th generation template brain at 2 mm resolution FIGURE 1 | The top panel shows the study protocol. The middle panel shows the display and timings for the Breath Hold (BH) task, Cued Deep Breathing (CDB) task, and REST sections. Visual instructions for each task were displayed on a monitor during scanning. For the BH task, IN and OUT instructions alternated for 3 s each, with a countdown from 24 s. Subjects ended on an exhale before holding and were instructed to do another exhale after holding. "Recover" is a period of free breathing. For the CDB task, IN and OUT instructions alternated for 2 s and subjects were told to take fast, deep breaths. The bottom panel displays the P ET CO 2 hrf traces (mmHg change from baseline) and GM-BOLD traces (% change from mean) for each fMRI acquisition. Each acquisition included 8 min of REST (gray), one acquisition preceded this with a BH task (blue) and one acquisition with a CDB task (green). From these three scans, five data segments of equal length were created. The first 10 volumes at the start of each data segment are not used (discarded to allow the signal to achieve a steady state of magnetization for BH + REST, CDB + REST and REST segments, and therefore matched for REST BH and REST CDB ) resulting in 390 volumes for each data segment. Thick lines represent group means, and thin lines represent each subject. (Grabner et al., 2006). These hemisphere masks were used alongside the Harvard Oxford Cortical Atlas (maxprob-thr25-2mm; Desikan et al., 2006), which parcellates the cortex into 48 regions, to make an atlas with 96 cortical parcels (HarvOx_96). The fsl_anat function outputs a linear transformation matrix from T1 space to MNI space; this was inverted and the HarvOx_96 atlas file was linearly transformed to T1 space for each subject using FSL FLIRT (Jenkinson and Smith, 2001;Jenkinson et al., 2002), with nearest neighbor interpolation and 12 degrees of freedom. The SBRef volume from the middle (second) fMRI dataset was registered to the preprocessed T1weighted image and all CVR maps were linearly transformed from fMRI space to T1-space for each subject. This resulted in the HarvOx_96 atlas file, CVR maps (No-Opt and Lag-Opt) and CBF maps all in T1-weighted space for each subject (Figure 2).

Between-Subject Correlation of Gray Matter Average Baseline Cerebral Blood Flow and Blood Oxygenation Level Dependent Cerebrovascular Reactivity
A representative GM CVR value and GM CBF value was computed by taking a median value over voxels inside the GM mask. For each of the five data segments, and for each optimization scheme (No-Opt and Lag-Opt), a Pearson's correlation analysis was performed between GM bCBF and BOLD-CVR values. Correlation plots and statistical outputs were created with the R packages ggplot2 (Wickham, 2016) and ggpubr FIGURE 2 | Example data files for one subject, in T1-weighted space. Top panel shows the Harvard-Oxford Cortical atlas: the initial 48 parcels were split into left and right hemisphere parcels making a total of 96. The middle panel shows CVR maps, not optimized for lag (No-Opt) and optimized for lag (Lag-Opt), for the BH + REST data segment. The bottom panel shows baseline cerebral blood flow (bCBF) maps. (Kassambara, 2020) with the "ggscatter" function. The Shapiro-Wilk test was used to ensure normality of variables.

Within-Subject Spatial Correlation of Baseline Cerebral Blood Flow and Blood Oxygenation Level Dependent Cerebrovascular Reactivity Across Cortical Regions
First, an average value of BOLD-CVR and bCBF within each atlas parcel (HarvOx_96, see Section "Transformation of Atlas and Cerebrovascular Reactivity Parameter Files to T1-Weighted Space") was computed. Second, for each analysis combination and each subject, a spatial correlation (Spearman's rank) between bCBF and BOLD-CVR was computed using custom MATLAB code (MathWorks, R2018b). The effect of data segment and lag optimization on bCBF-CVR spatial correlation values was tested with a repeated measures ANOVA, run with the R package permuco (Frossard and Renaud, 2021) with the "aovperm" function. Null distributions were created via 100,000 permutations of the original data, which therefore do not depend on Gaussian and sphericity assumptions. When investigating simple main effects ("emmeans" package) and performing multiple comparisons, p-values were false discovery rate (FDR) corrected, and then compared against an alpha of 0.05 to determine significance.

Within-Subject Correlation of Resting-State Metrics and Baseline Cerebral Blood Flow Across Cortical Regions
Results from our previous work (Stickland et al., 2021) and from the analyses described in Sections "Between-Subject Correlation of Gray Matter Average Baseline Cerebral Blood Flow and Blood Oxygenation Level Dependent Cerebrovascular Reactivity" and " Within-Subject Spatial Correlation of Baseline Cerebral Blood Flow and Blood Oxygenation Level Dependent Cerebrovascular Reactivity Across Cortical Regions" show the difficulty of estimating BOLD-CVR with resting-state data, within a lagged-GLM framework using a P ET CO 2 regressor. We hypothesized this to be signal-to-noise driven: if intrinsic low-frequency oscillations, associated with neural activity or other physiological sources, are similar or greater magnitude to the low frequency fluctuations induced by P ET CO 2 , this could lead to poor coupling between P ET CO 2 and fMRI signal changes. However, literature suggests that metrics derived directly from BOLD resting-state fluctuations may capture CVR effects (Chen and Gauthier, 2021;Pinto et al., 2021), albeit not in quantitative units, therefore making comparisons between subjects or between time-points more limited. We performed a simple exploratory analysis looking at the spatial relationship between resting-state metrics and bCBF. Minimally pre-processed BOLD-fMRI data (brain extraction, volume registration) from the three resting-state data segments (REST, REST BH and REST CDB ) were input to AFNI's 3dRSFC function (including quadratic detrending of input data, bandpass filtering between 0.01 and 0.1 Hz and spatial smoothing with a 4 mm full width half maximum (FWHM) kernel). This function computed maps of resting-state fluctuation amplitude normalized to the mean (mRSFA), amplitude of low-frequency fluctuation (ALFF; Zang et al., 2007) and fractional ALFF (fALFF; Zou et al., 2008). ALFF and RSFA (Kannurpatti and Biswal, 2008) are similar metrics; the normalized version of RSFA (mRSFA) better distinguishes this metric from ALFF (Liu et al., 2013b;Golestani et al., 2016). The same spatial correlations with bCBF were performed for each subject as explained in Section "Within-Subject Spatial Correlation of Baseline Cerebral Blood Flow and Blood Oxygenation Level Dependent Cerebrovascular Reactivity Across Cortical Regions". Figure 3 shows the relationship between GM bCBF and GM BOLD-CVR across subjects, for each data segment and lag optimization scheme. After lag optimization (and with influential points removed), all correlations between bCBF and CVR are positive: a higher bCBF tends to co-occur with a higher CVR across subjects. Importantly, the only significant correlations with bCBF were found with lag optimized CVR values derived from data segments with breathing tasks included (i.e., BH + REST and CBD + REST). However, no correlations were found to be statistically significant after FDR correction across the 10 tests. These results are consistent across different GM masking options (Supplementary Table 1). Figure 4 illustrates the spatial correlation, for each subject, across the 96 cortical parcels of the atlas. Single subject correlations that are significant (outside the pink shaded box in panel B) are mostly positive. Panel C of Figure 4 shows how bCBF and BOLD-CVR correlations changed due to lag optimization. This change was most consistent for BH + REST data, with 8/9 subjects showing an increased (positive) correlation after lag optimization. For the other segments, changes with lag optimization are more variable, and there are more negative correlations. The repeated measures ANOVA showed no significant interaction effect (F(4,32) = 0.77, p = 0.5552) and two significant main effects: correlation values were significantly higher for Lag-Opt versus No-Opt (F(1,8) = 6.09, p = 0.0385) and there was a significant effect of data segment on correlation values (F(4,32) = 3.84, p = 0.0117). BH + REST showed significantly higher values than REST CDB (see Supplementary Table 2 for all pairwise comparisons).

Within-Subject Spatial Correlation of Baseline Cerebral Blood Flow and Blood Oxygenation Level Dependent Cerebrovascular Reactivity Across Cortical Regions
Two extreme outlier points were identified (greater than 3 × IQR below the first quartile or above the third), indicated by an asterisk in Figure 4. When removing these two subjects the same results were found: no significant interaction effect (F(2,24) = 0.30, p = 0.8732), and significant main effects of optimization scheme (F(1,6) = 8.34, p = 0.02682) and data segment (F(2,24) = 10.83, p = 0.00006). When replacing these outlier values with the mean value for that cell of the design, the same results were also found: no significant interaction (F(4,32) = 0.46, p = 0.7661), and significant main effects of optimization scheme (F(1,8) = 10.67, p = 0.0111) and data segment (F(4,32) = 7.40, p = 0.00031). Multiple ANOVAs were explored due to the challenges with removing outliers, particularly for within-subject designs. Figure 5 illustrates that in the REST data segment, all single subject spatial correlations between bCBF and ALFF, fALFF and mRSFA were positive and significant at the single subject level, with fALFF showing slightly higher correlation strength with bCBF compared to ALFF and mRSFA. ALFF and mRSFA showed extremely similar maps, and therefore also similar correlation values with bCBF, which is logical because a correlation analysis is not sensitive to differences in absolute scale. A similar pattern of results was also seen for the REST BH and REST CDB data segments (Supplementary Figure 1). It is important to note that we generated outputs with and without smoothing (Supplementary Figure 1); the correlations are notably different for ALFF and mRSFA, showing much lower correlation values without smoothing, which are no longer significant at the singlesubject level for the majority of subjects. Smoothing is a typical processing step in resting-state fMRI analyses, yet the appropriate kernel size is hard to determine. Using a Gaussian kernel size of two or three times the voxel size is a common recommendation (Alahmadi, 2021;Cumming et al., 2012;Liu et al., 2017;Pajula and Tohka, 2014); thus, we smoothed with FWHM size of 4 mm, double our voxel size.

DISCUSSION
We assessed the impact of methodological factors on the correlation between bCBF and BOLD-CVR. We present complimentary results to those in our previous manuscript (Stickland et al., 2021), suggesting that adding a simple breathing FIGURE 3 | Between-subject Pearson correlations between GM CBF values and GM BOLD-CVR values, shown for each data segment (1-5) and for CVR values with no lag optimization (No-Opt, black circles) and CVR values with lag optimization (Opt, pink triangles). The CBF and CVR values are medians over a GM mask and each dot represents one participant. On the right, statistical tests are shown with and without influential points removed. Influential points were classed as points with a Cook's distance over 4/n, with n being the number of subjects. P-values are not corrected for multiple comparisons. The gray shaded regions around the fit lines indicate the 95% confidence interval of the correlation coefficient (fit line and confidence interval estimation did not include influential points, but influential points are still plotted alongside).
task to a resting-state BOLD fMRI dataset, alongside lagoptimization of a P ET CO 2 regressor within CVR modeling, improves the relationship between BOLD-CVR and bCBF. We also report exploratory results of positive spatial correlations between bCBF and ALFF, mRSFA and fALFF, for each of the three resting-state segments. We first discuss the impact of these methodological factors in more detail (Section "The Impact of Methodological Factors on the Correlation Between Baseline Cerebral Blood Flow and Blood Oxygenation Level Dependent Cerebrovascular Reactivity"), followed by exploring reasons for the observed positive relationship between bCBF and BOLD-CVR and the resting-state metrics (Section "A Positive Correlation Between Baseline Cerebral Blood Flow and Blood Oxygenation Level Dependent Cerebrovascular Reactivity"), and then consider the relevance and implications for fMRI research (Section " Relevance and Implications for Functional Magnetic Resonance Imaging Research").

The Impact of Methodological Factors on the Correlation Between Baseline Cerebral Blood Flow and Blood Oxygenation Level Dependent Cerebrovascular Reactivity
For variability between-subjects, a significant positive correlation between an average GM BOLD-CVR value and an average GM bCBF value was only found for BOLD-CVR data that included a breathing task and when analysis included optimization for hemodynamic lag effects. With influential points removed, the strength of the correlation always increased after lagoptimization, for all comparisons. For within-subject spatial correlations between bCBF and BOLD-CVR in cortical parcels, the bCBF and BOLD-CVR relationship was mostly positive (the negative correlations were only seen in resting-state data with no breathing tasks or in results without lag-optimization). At the group level, spatial correlations were significantly higher after lag-optimization of the CVR values. This suggests that these two methodological choices, the inclusion of breathing data and the implementation of lag-optimization, result in more physiologically accurate CVR estimates.
We have shown previously (Stickland et al., 2021) that adding breathing tasks to resting state fMRI paradigms makes CVR mapping more robust: if BOLD signal changes related to P ET CO 2 are small in amplitude, as may occur if natural fluctuations in P ET CO 2 are small during resting-state paradigms, it can be hard to distinguish them from other physiological, artifactual, or neuronally-driven fluctuations that occur at the same low frequencies. Our observations in this work suggest that BOLD-CVR estimates derived from resting-state data also show a much more variable relationship with bCBF. The timing of the blood flow response to a change in the pressure of arterial CO 2 (as measured here with end-tidal CO 2 recordings) will not be the same for each equipment setup, each participant, and each area of the brain due to methodological and physiological factors. As well as the data type to be modeled (resting state only vs. added breathing task), accounting for voxel-wise hemodynamic timings when estimating BOLD-CVR seems to lead to more physiologically grounded CVR values.
For segments that contained breathing tasks, BH + REST showed more consistent patterns than CDB + REST. Specifically, for the between-subject correlations of GM average values, both BH + REST and CDB + REST showed positive correlations FIGURE 4 | The association between bCBF and CVR values in 96 cortical parcels was computed for each subject and analysis combination (panel A). A Fisher's Z transformation is applied to the correlation values to allow group statistical testing; panel (B) visualizes these transformed correlations for the five segments (columns), both for CVR values with no lag optimization (No-Opt) and CVR values with lag optimization (Opt). The shaded pink box in panel B represents correlations that were not significant at the single-subject level, based on 96 datapoints and a critical value of 1.96 (p < 0.05, two-tailed). When generating these correlation values for each subject, outlier parcels (based on Cook's distance) were not included. The number of parcels removed from each correlation analysis averaged 5.3 out of 96 (see Supplementary Table 3). As removing parcels reduces the degrees of freedom, the dotted pink line shows the adjusted critical value (maximum across participants and data segments shown for reference). Panel (C) shows the change in absolute correlation strength due to lag optimization of the CVR values.
between GM bCBF and BOLD-CVR that increased after lag-optimization, but there was one subject outlier for the CDB + REST correlation analysis. For within-subject spatial correlations, the BH + REST segment had higher correlation values on average compared to CDB + REST (though this was not a significant pairwise difference), and a more consistent change with lag optimization: 8/9 subjects showed an increased (positive) correlation after lag optimization for BH + REST, whereas this was only the case for 7/9 subjects in the CDB + REST segment, with one of these seven subjects having nearly zero correlation after lag-optimization. It is unclear what drives these differences; it simply may be the case that we can more consistently get accurate CVR estimates with the BH + REST compared to the CDB + REST segment (Stickland et al., 2021), and therefore more FIGURE 5 | The association between bCBF and each resting-state metric, in 96 cortical parcels, was computed for each subject, using data from the REST data segment. Maps from one example subject are displayed. A Fisher's Z transformation is applied to the correlation values before averaging across subjects. The shaded pink box represents correlation values that would not be significant at the single-subject level, based on 96 datapoints and a critical value of 1.96 (p < 0.05, two-tailed). When generating these correlation values for each subject, outlier parcels (based on Cook's distance) were not included. As removing parcels reduces the degrees of freedom, the dotted pink line shows the adjusted critical value (maximum across participants and resting-state metric shown for reference).
consistent correlation with bCBF can be observed. The reasons for this cannot be fully addressed with our small sample and this study was not designed for an effective comparison between these two breathing tasks, yet some interpretations can be explored. Differences in CVR mapping could be due to the BH task being slightly longer, including three cycles of hypercapnia versus two cycles of hypocapnia for CDB, as well as the inclusion of pacedbreathing sections as opposed to free breathing. Alternatively, the capacity for vasodilation (driven by BH) may be physiologically coupled to bCBF differently compared to the capacity for vasoconstriction (driven by a CBD task); previous work has demonstrated that CVR can be modulated differently by baseline vascular tension (Bright et al., 2011;Halani et al., 2015) and that vasoconstrictive and vasodilatory responses can be differentially affected in pathology (Zhao et al., 2009).
It is not clear why the REST CDB segment had significantly lower spatial correlations between bCBF and BOLD-CVR than the other two REST segments, irrespective of lag-optimization. It could possibly be due to residual vascular after-effects following the hypocapnic breathing task stimulus (Bright et al., 2011). The CDB task had fewer visual instructions than the BH task, so if neurovascular after-effects played a role, we would expect this to be more relevant for the REST BH data segment as opposed to REST CDB . Furthermore, the REST CDB data segment that was analyzed was preceded by 55 seconds of the subject breathing normally (the last 43 seconds of the "breathe normally" task cue, see Figure 1, and the first 12 seconds of the REST CDB section not analyzed during CVR modeling). We previously found no clear evidence for obvious differences in common restingstate metrics (RSFA, ALFF, fALFF, LFCD) between all three resting-state data segments (Stickland et al., 2021), though further investigation of neural and vascular after-effects of breathing tasks may be warranted.

A Positive Correlation Between Baseline Cerebral Blood Flow and Blood Oxygenation Level Dependent Cerebrovascular Reactivity
There are many studies in the literature including both bCBF and BOLD-CVR values for healthy young populations, but few studies that directly assess the covariance of these values. In conflict with our results, a PET study of healthy subjects, of similar age to this study, showed a significant negative correlation between bCBF and CBF-CVR (with hypocapnia and hypercapnia) in multiple brain regions (Ito et al., 2008). This difference could be driven by the fact they assessed CBF-CVR, as opposed to BOLD-CVR, which we discuss further in a subsequent section. A recent study  with a sample of similar age to ours investigated sex differences in bCBF and BOLD-CVR, showing that female participants had higher bCBF (consistent with what we see in Supplementary  Figure 2) but no difference in BOLD-CVR. However, they also reported a negative correlation between bCBF and BOLD-CVR, using average values from four different brain lobes (frontal, temporal, parietal, and occipital). They concluded that low and high autoregulation efficiencies exist at high and low bCBF values, respectively. This is again contrary to our results, where we see positive correlations between these metrics. There are key methodological and analytical differences between this study and ours: average values were derived from 4 brain lobes compared to the 96 cortical parcels in this study; they did not record P ET CO 2 values, precluding characterization of BOLD-CVR in units of%BOLD/mmHg; the BOLD response to the breath-hold was modeled with the convolved task design and hemodynamic timings were not corrected on a regional level; CBF and BOLD fMRI data were smoothed with an 8 mm FWHM Gaussian kernel whereas we did not apply any spatial smoothing except the adaptive spatial smoothing applied with FSL BASIL's toolbox. Finally, the focus of their study was to compare sex differences in bCBF and BOLD-CVR; their CBF and CVR correlations are displayed as group averages with large standard deviations, and it is therefore hard to distinguish the contribution of individual variability to this main result.

Aging Effects
Previous publications reporting or commenting on a positive correlation between bCBF and BOLD-CVR across subjects often explain this in the context of cerebrovascular aging effects: in elderly adults, a smaller capacity for vasodilation is often seen alongside lower bCBF (Marstrand et al., 2002;Lu et al., 2011;Leoni et al., 2017;Tsvetanov et al., 2020). Developmental age is an important factor to consider; one study examined how BOLD-CVR and bCBF correlated with age in a sample ranging from 9 to 30 years (Leung et al., 2016). BOLD-CVR and bCBF only changed in the same way with age after 14.7 years (both decreased with age), whereas between 9 and 14.7 years bCBF decreased with age and CVR increased with age, resulting in a negative relationship between these parameters. We cannot rule out age effects as a potential causal explanation for the bCBF-CVR relationships in our data, and we do see negative trends with age and bCBF and negative trends with age and BOLD-CVR that agree with the literature (Supplementary Figure 2). However, it is challenging to interpret our findings within the framing of aging effects or cerebrovascular impairment given our small sample of healthy younger adults with a small age range (23-35 years), and this was not the focus of the study. Therefore, alternative interpretations separate from aging effects are explored next.

Neurovascular and Metabolic Individual Differences
Under the assumption that our vasoactive breathing task stimuli are isometabolic, the BOLD-CVR response should reflect the CBF-CVR response (Zhou et al., 2015). However, there are many physiological factors that can contribute to a BOLD signal change other than a CBF change (Kim and Ogawa, 2012), and the assumption of an isometabolic state under hypercapnia may not always be valid (Driver et al., 2017;Yablonskiy, 2011;Deckers et al., 2021;Xu et al., 2011). Furthermore, BOLD-CVR and CBF-CVR may have differing dependencies on the baseline condition (Stefanovic et al., 2006;Halani et al., 2015) and on the type of vascular stimulus (Halani et al., 2015). It is therefore possible that individual variability in BOLD-CBF coupling during breathing tasks could influence our results.
Healthy individual differences in the functioning of the neurovascular unit (Schaeffer and Iadecola, 2021) could also affect both bCBF and BOLD-CVR in similar ways, leading to their positive coupling. To contextualize our results within commonly derived resting-state BOLD metrics, separate from our lagged-GLM framework, we performed an exploratory analysis to assess the relationship between resting-state metrics and bCBF. We found weak to moderate positive spatial correlations with bCBF for both ALFF and mRSFA (these are directly proportional measures). The contribution of neural or vascular factors driving voxel-wise ALFF/mRSFA values is not clear, and they have been interpreted and applied as both neural and vascular in the literature (Chen and Gauthier, 2021). Interestingly, fALFF had a stronger positive correlation with bCBF in our data compared to ALFF/mRSFA (regardless of smoothing, Supplementary  Figure 1), and fALFF is thought to be more neuronally specific than ALFF. A recent paper discusses and investigates the complexities of these physiological interpretations (Deng et al., 2022), and their results suggest ALFF is more influenced by venous vascular and microvascular density (Vigneau-Roy et al., 2014), and fALFF is reflective of underlying metabolic demand. Therefore, the positive spatial correlation between fALFF and bCBF could reflect more metabolically active regions at rest requiring more blood flow. The degree of coupling between resting-state BOLD and CBF signals also relates to the macrovascular volume fraction (i.e., BOLD and CBF are less coupled in voxels near larger vessels; Tak et al., 2014). It is therefore possible that the positive coupling of bCBF with BOLD-CVR, ALFF/mRSFA and fALFF could be caused by overlapping, yet distinct physiological drivers (i.e., BOLD driven effects, CVR driven, metabolism driven).

Non-linear Relationship Between Cerebral Blood Flow and Partial Pressure of End Tidal Carbon Dioxide
This positive bCBF and BOLD-CVR coupling could reflect the non-linear relationship between CBF and P ET CO 2 : the sigmoidal dose-response curve combined with lowered baseline P ET CO 2 levels associated with paced breathing between breath-holds, or due to sustained effects of hyperventilation in the CDB task, may strongly influence the magnitude of CBF change during breathing tasks [for example, see Tancredi and Hoge (2013)]. However, there is evidence that CVR studies can still assume to operate in the linear portion of the dose-response curve if baseline P ET CO 2 values are kept within 30 -45 mmHg (Tancredi and Hoge, 2013) as is the case with our data (Supplementary Table 4, Part B). Nevertheless, CVR studies using breathing tasks should take care to avoid causing undue hypocapnia during baseline periods, whilst attempting to achieve stability with paced breathing, as this could bias CVR measurements.
Finally, a recent study aimed to identify physiological factors that contribute to the large observed between-subject and between-session variability seen within BOLD-CVR data (Hou et al., 2020); the authors reported a significant negative association between BOLD-CVR and baseline P ET CO 2 , across subjects. The relationship between baseline P ET CO 2 and CVR is likely mediated by bCBF, yet when performing exploratory correlation analyses, we did not find significant evidence that bCBF or BOLD-CVR was related to baseline P ET CO 2 in our data (Supplementary Table 4). We interpret these results cautiously considering our much smaller sample size compared to the previous work (Hou et al., 2020). Furthermore, our use of a nasal cannula is valid to sample P ET CO 2 (Kim et al., 1997;Yosefy et al., 2004;Nayak et al., 2019), particularly in a healthy sample, but when estimating absolute baseline P ET CO 2 for any individual this is likely to be underestimated due to the contribution of room air during sampling.

Relevance and Implications for Functional Magnetic Resonance Imaging Research
A key motivation behind understanding the relationship between bCBF and BOLD-CVR is to facilitate appropriate normalization analysis strategies, as has been done for task-based BOLD fMRI. For example, previous work has corrected task-induced BOLD signals by covarying out bCBF (Krishnamurthy et al., 2020) or CVR (Cohen et al., 2004;Liu et al., 2013a) to account for varying neurovascular coupling or vascular physiology. However, based on the conflicting evidence presented in existing literature and the results we report here, we suggest that the coupling between bCBF and BOLD-CVR needs to be understood better before such normalization strategies become common-place, and more consensus needs to exist in the healthy population literature. Though we report a significant positive correlation between bCBF and BOLD-CVR, they only share a proportion of their variance. For example, GM bCBF and lag-optimized GM BOLD-CVR shared 53% and 61-67% of between-subject variance for BH + REST and CDB + REST data, respectively. For the within-subject spatial analysis, the average variance shared between bCBF and BOLD-CVR was 42% and 33% for BH + REST and CDB + REST data, respectively. Therefore, it will be important to consider both these metrics when studying the cerebrovasculature in healthy populations, whilst modeling their effects in ways that account for shared variance. The validity of fMRI normalization strategies or covariate analyses will depend on how physiologically accurate the estimates of bCBF and BOLD-CVR are. Many factors will affect the accuracy; in this paper and our previous work (Stickland et al., 2021) we have shown evidence for two factors affecting CVR estimates (lag-optimization and using breathing task data vs. restingstate data only).
Studying the relationship between bCBF and BOLD/CBF-CVR becomes more interesting, yet more complicated, in situations where they may become uncoupled. In populations where we know or suspect bCBF or CVR to be altered, causing an uncoupling of neuronal activation and blood flow (Para et al., 2017), it is even more important to understand the covariance between bCBF and CVR. However, normalizing fMRI results using these metrics in situations where they are uncoupled or coupled in atypical ways is problematic, particularly if causality is not clear (i.e., whether an impairment of bCBF drives an impairment of CVR, or vice versa). For example, in children with sickle cell disease there are reports of elevated bCBF but lower BOLD and CBF CVR (Kosinski et al., 2017;Václavů et al., 2019); in cognitively normal older adults, greater aortic stiffness has been related to lower regional bCBF and higher CBF-CVR (Jefferson et al., 2018); in healthy adults with a normal fitness range, there are reports of aerobic fitness being associated with lower bCBF but greater CBF-CVR (Foster et al., 2020). There are many more examples in the literature where bCBF and CVR do not show the positive coupling observed in this study. Ultimately, understanding the relationship between bCBF and BOLD/CBF-CVR estimates in healthy populations, and what factors affect this relationship, will help us better assess the extent to which they are uncoupled in pathology, and better interpret the nature of this uncoupling.

Study Limitations
This study is limited in its conclusions by its small sample size. We chose to visualize each subject's data points to be transparent regarding the full data variance, rather than only reporting group summary measures. We only used a single post label delay for our pCASL acquisition, which did not allow us to estimate and account for arterial arrival times in bCBF quantification. However, the post-labeling delay in our acquisition was chosen to be long enough for labeled spins to reach the imaged tissue in healthy young brains (Alsop et al., 2015), and therefore perfusion quantification should be accurate. For the vasoactive stimuli, we modulated P ET CO 2 values with a breathing task or characterized resting-state fluctuations in P ET CO 2 , as opposed to delivering gas inhalation challenges. Such gas challenges, particularly when administered by feed-forward gas blending systems (Slessarev et al., 2007), would allow more precise and repeatable P ET CO 2 targeting while keeping P ET O 2 values more stable. This methodology could be implemented in future work to improve mapping of BOLD-CVR and further probe its relationship with bCBF. A limitation already discussed in the section "Neurovascular and Metabolic Individual Differences" is that we are using the BOLD signal contrast to reflect CBF changes. Future work should consider the measurement of CBF-CVR as well as BOLD-CVR, when probing the bCBF-CVR relationship across individuals and brain regions. Although it is possible to measure CBF-CVR with breathing task designs (Solis- Barquero et al., 2021), controlled gas inhalation challenges will greatly facilitate accurate CBF-CVR measurements through repeated arterial spin labeling acquisitions that demand sustained stable P ET CO 2 modulations.

CONCLUSION
In summary, we correlated two separately acquired metrics of cerebrovascular health, bCBF and BOLD-CVR (with BOLD-CVR estimates from 3 separate fMRI acquisitions), in 9 healthy individuals. These metrics showed stronger and more consistent positive correlations when CVR was modeled with breathing task data as opposed to restingstate data only, and when CVR timing was accounted for. Using a vasodilatory breathing-task (BH) produced more consistent results than a vasoconstrictive breathing-task (CDB). The stronger relationship between bCBF and BOLD-CVR suggests more physiologically valid CVR estimates when such methodological elements are incorporated. Future research should investigate the basis of the observed positive relationship between bCBF and BOLD-CVR in larger normative samples and across the lifespan. Understanding how baseline vascular physiology relates to dynamic vascular processes is particularly important in clinical cohorts that present with altered vascular and/or metabolic baselines or impaired cerebrovascular reserve.

DATA AVAILABILITY STATEMENT
The datasets presented in this article are not readily available because of restrictions of the ethical approval that they were collected under. However, analysis derivatives that are not included in this manuscript may be provided, on request, within ethical guidelines. Requests to access the datasets should be directed to MB, molly.bright@northwestern.edu.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Northwestern University's Institutional Review Board. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
RS: conceptualization, methodology, software, formal analysis, investigation, data curation, writing (OD), writing (RE), visualization, and project administration. KZ: methodology, investigation, and writing (RE). SM and CC-G: methodology and writing (RE). MB: conceptualization, methodology, software, investigation, resources, writing (RE), supervision, project administration, and funding acquisition. All authors contributed to the article and approved the submitted version.