Cerebrovascular reactivity assessment with O2-CO2 exchange ratio under brief breath hold challenge

Background Hypercapnia during breath holding is believed to be the dominant driver behind the modulation of cerebral blood flow (CBF). However, increasing evidence show that mild hypoxia and mild hypercapnia in breath hold (BH) could work synergistically to enhance CBF response. We hypothesized that breath-by-breath O2-CO2 exchange ratio (bER), defined as the ratio of the change in partial pressure of oxygen (ΔPO2) to that of carbon dioxide (ΔPCO2) between end inspiration and end expiration, would be able to better correlate with the global and regional cerebral hemodynamic responses (CHR) to BH challenge. We aimed to investigate whether bER is a more useful index than end-tidal PCO2 to characterize cerebrovascular reactivity (CVR) under BH challenge. Methods We used transcranial Doppler ultrasound (TCD) to evaluate CHR under BH challenge by measuring cerebral blood flow velocity (CBFv) in the middle cerebral arteries. Regional changes in CHR to BH and exogenous CO2 challenges were mapped with blood oxygenation level dependent (BOLD) signal changes using functional magnetic resonance imaging (fMRI). We correlated respiratory gas exchange (RGE) metrics (bER, ΔPO2, ΔPCO2, end-tidal PCO2 and PO2, and time of breaths) with CHR (CBFv and BOLD) to BH challenge. Temporal features and frequency characteristics of RGE metrics and their coherence with CHR were examined. Results CHR to brief BH epochs and free breathing were coupled with both ΔPO2 and ΔPCO2. We found that bER was superior to either ΔPO2 or ΔPCO2 alone in coupling with the changes of CBFv and BOLD signals under breath hold challenge. The regional CVR results derived by regressing BOLD signal changes on bER under BH challenge resembled those derived by regressing BOLD signal changes on end-tidal PCO2 under exogenous CO2 challenge. Conclusion Our findings provide a novel insight on the potential of using bER to better quantify CVR changes under BH challenge.


Background
Hypercapnia during breath holding is believed to be the dominant driver behind the modulation of cerebral blood flow (CBF). However, increasing evidence show that mild hypoxia and mild hypercapnia in breath hold (BH) could work synergistically to enhance CBF response. We hypothesized that breath-by-breath O 2 -CO 2 exchange ratio (bER), defined as the ratio of the change in partial pressure of oxygen (ΔPO 2 ) to that of carbon dioxide (ΔPCO 2 ) between end inspiration and end expiration, would be able to better correlate with the global and regional cerebral hemodynamic responses (CHR) to BH challenge. We aimed to investigate whether bER is a more useful index than end-tidal PCO 2 to characterize cerebrovascular reactivity (CVR) under BH challenge.

Methods
We used transcranial Doppler ultrasound (TCD) to evaluate CHR under BH challenge by measuring cerebral blood flow velocity (CBFv) in the middle cerebral arteries. Regional changes in CHR to BH and exogenous CO 2 challenges were mapped with blood oxygenation level dependent (BOLD) signal changes using functional magnetic resonance imaging (fMRI). We correlated respiratory gas exchange (RGE) metrics (bER, ΔPO 2 , ΔPCO 2 , endtidal PCO 2 and PO 2 , and time of breaths) with CHR (CBFv and BOLD) to BH challenge. Temporal features and frequency characteristics of RGE metrics and their coherence with CHR were examined.

Results
CHR to brief BH epochs and free breathing were coupled with both ΔPO 2 and ΔPCO 2 . We found that bER was superior to either ΔPO 2 or ΔPCO 2 alone in coupling with the changes of

Introduction
Breath hold challenge has been used in the clinical setting as a simple vasoactive stimulus for the assessment of cerebrovascular reactivity (CVR) [1,2] in patients with carotid artery diseases [3][4][5] as well as brain tumors [6]. Its use for CVR assessment with transcranial Doppler sonography (TCD) measurement was first demonstrated by Ratnatunga and Adiseshiah [3].
The cerebral blood flow velocities (CBFv) were often measured in the middle cerebral arteries which supply most parts of the brain. The time of breaths (ToB) had often been taken as an indicator of the strength of the vasoactive stimulus to induce changes of cerebral blood flow (CBF). The ratio of CBF to ToB was known as a breath hold index to evaluate CVR [1,2,4,5,7]. Although TCD offers high temporal resolution to evaluate cerebrovascular responses without the concern of aliasing high frequency hemodynamic signal into the low frequency range, it does not provide regional information. Therefore regional CVR mapping with blood oxygen level-dependent (BOLD) signal changes measured by functional magnetic resonance imaging (fMRI) was used. BOLD-fMRI was used instead of arterial spin labeling (ASL) in MRI perfusion studies due to the low contrast to noise ratio and the low temporal resolution of the ASL technique [8]. ASL image acquisition at a temporal resolution of 4 seconds will under-sample the brain responses within respiratory cycle of 4-6 seconds. Hypercapnia or increased end-tidal partial pressure of CO 2 (P ET CO 2 ) was commonly measured as a surrogate for the increased arterial CO 2 level to evaluate CVR [9][10][11]. Hypoxia or changes of end-tidal partial pressure of O 2 (P ET O 2 ) were seldom used to account for CBF changes during breath holding partly due to the common belief that the vasodilatory effect of increased arterial partial pressure of CO 2 (PCO 2 ) dominates that of decreased arterial partial pressure of O 2 (PO 2 ). Such a belief stemmed from high altitude studies and lab-controlled low oxygen environment reporting that significant CBF changes only happened in hypoxia with arterial PO 2 (PaO 2 ) going down to approximately 50 mmHg [12][13][14]. It should be noted that reports of large CBF increase in response to strong hypoxia were often accompanied by hypocapnia due to hyperventilation in the studies of either high altitudes or lab-controlled low oxygen environment (~8-13% O 2 ) [12][13][14]. The case is different with breath hold because hypoxia is accompanied by hypercapnia and not by hypocapnia during breath holding.
The role of both hypoxia and hypercapnia in breath hold deserves to be examined because arterial PO 2 and PCO 2 had been reported to work synergistically at normoxia to interact with peripheral chemoreceptors [15][16][17][18][19]. In spontaneous breathing, blood gas levels of O 2 and CO 2 are optimized by the feedback control of ventilation via chemoreflexes [20] to regulate blood flow and oxygen delivery to the brain as part of a vital homeostatic process [17]. The feedback loops include interaction between central chemoreceptors at the brainstem [21,22] and peripheral chemoreceptors at the carotid body [23,24]. While most of the studies on the role of PO 2 to stimulate peripheral chemoreceptors at the carotid body had been focused on strong hypoxia [17] where chemoreceptor activities rose quickly in a hyperbolic fashion, Biscoe et al. [15] showed that peripheral chemoreceptor activities could be observed from normoxia to hyperoxia up to arterial PO 2 level of 190mmHg and beyond. Lahiri et al. [18] reported that the stimulus thresholds of arterial PO 2 and PCO 2 for peripheral chemoreceptors were largely interdependent under the normoxic condition where a drop in arterial PO 2 was routinely accompanied by increased chemoreceptor activities as well as an enhanced sensitivity of carotid chemoreceptors to arterial PCO 2 . Several authors [15][16][17][18][19] reported similar findings for normoxic as well as hypoxic conditions. Ventilation, like chemoreceptors, also becomes more sensitive to PCO 2 with a slight decrease in PO 2 within the range of normoxia (90-110 mmHg). The same interaction between change in PO 2 and change in peripheral chemoreceptor activities during spontaneous breathing is also assumed to take place during breath hold.
In terms of the relationship between modulated chemoreceptor activities and cerebral hemodynamic responses, previous studies reported that apnea-induced hypoxia and hypercapnia caused chemoreceptor-mediated central vasodilation and concurrent peripheral vasoconstriction to conserve oxygen delivery to the brain [25], leading to an increase in cerebral blood flow (CBF) and a decrease in peripheral oxygen saturation [26]. Enhanced hemodynamic responses to mild hypoxia with maintained eucapnia or in the presence of hypercapnia in humans had also been demonstrated in the work by Shapiro et al. [27] and Mardimae et al. [28]. Those authors reported progressive CBF increase in response to small steps of serial reductions of PO 2 starting from normoxia with maintained eucapnia [27] as well as with slight hypercapnia of around 45 mmHg [28]. The findings in these studies suggested that mild hypercapnia could increase the sensitivity of the CBF response to a very mild level of hypoxia and the ranges of mild PO 2 and PCO 2 changes reported are achievable by breath hold.
Beyond the belief of CO 2 being the dominant effect to account for CBF changes, the perception that effect of O 2 on CBF simply reflects the effect of CO 2 is often associated with another belief that changes of PO 2 and PCO 2 in respiration are mirror image of each other [29]. However, respiratory data acquired on 12 subjects in the study by Lenfant [30] showed that the time courses of O 2 and CO 2 were interdependent but not redundant in terms of their temporal and frequency characteristics. Even without being redundant, such interdependence can still pose some problems for CVR analysis. There is a question of collinearity when interdependent changes of PCO 2 and PO 2 are included as predictor variables in the same regression model of CVR quantification. We are seeking a combination of PO 2 and PCO 2 that can properly characterize the interdependence between PO 2 and PCO 2 without creating a problem of collinearity in CVR analysis. The use of the breath-by breath O 2 -CO 2 exchange ratio described in the alveolar air equation [31,32] is a natural target to be investigated for breath hold challenge in the regression model.
We preferred the terms ΔPO 2 (inspired PO 2 -expired PO 2 ) and ΔPCO 2 (expired PCO 2inspired PCO 2 ) over the more commonly used P ET O 2 and P ET CO 2 because the breath-bybreath ΔPO 2 and ΔPCO 2 are the terms used in the alveolar air equation [31,32] to describe the change of gas partial pressure in systemic O 2 uptake and that in CO 2 release respectively. We also prefer to use ΔPO 2 /ΔPCO 2 as our breath-by-breath O 2 -CO 2 exchange ratio (bER). bER takes advantage of its ratio format to reduce the unwanted effects of ventilatory volume fluctuations due to isolated deep breaths common to ΔPO 2 and ΔPCO 2 measured. bER is mathematically equivalent to the reciprocal of the respiratory exchange ratio (RER) from alveolar air equation [31,32], except that bER is a dynamic breath-by-breath measurement and RER is a steady-state time-averaged measurement over a period of time. RER has been used to evaluate resting systemic metabolic rate [33][34][35]. Some technical differences do need to be mentioned between RER used in the literature and bER we used in this study. Traditionally, RER is derived by measuring the respiratory flow and the expired gases collected in Douglas bag connected to a closed circuit over several minutes. bER is derived here by measuring the inspired and expired gases with a nasal tubing at each breath.
In this study, we hypothesized that mild hypoxia and hypercapnia work synergistically to increase CBF under breath hold challenge. The primary objective was to evaluate our hypothesis that bER would be a more useful index than P ET CO 2 to correlate with the changes of CBFv and BOLD signals in the evaluation of cerebral hemodynamic responses to the breath hold challenge. To address the question of redundancy among respiratory gas exchange (RGE) metrics, we examined the correlations among bER, ΔPCO 2 , ΔPO 2 , P ET CO 2 , P ET O 2 and ToB in both TCD and MRI sessions. In the first part of the study, we correlated the time series of the RGE metrics (bER, ΔPCO 2 , ΔPO 2 and ToB) with those of CBFv and BOLD signal changes under the protocol of breath hold challenge. We then examined the temporal features and frequency characteristics of these RGE metrics as well as their coherence with CBFv and BOLD signal changes. In the second part of the study, we compared the regional CVR maps obtained by regressing BOLD signal changes on selected RGE metrics including bER, P ET CO 2 and ToB under breath hold challenge with regional CVR maps obtained by regressing BOLD signals changes on P ET CO 2 obtained under exogenous CO 2 hypercapnic challenge. The success of such association between bER and cerebral hemodynamic responses, in addition to offering a better physiological model to characterize cerebral hemodynamic responses under breath hold challenge, would also provide a novel insight in the study of brain-body interaction.

Participants
Seventeen volunteers aged from 22 to 48 years (mean = 31.18 years; SD = 8.78 years; 11M and 6F) were included. All of them were recruited by e-mail and poster placement within the Partners hospital network. They were screened to exclude neurological, mental and medical disorders and drug abuse. TCD and MRI scanning were performed in the Athinoula A. Martinos Center for Biomedical Imaging at the Massachusetts General Hospital of Partners HealthCare. All the experimental procedures were explained to the subjects, and signed informed consent was obtained prior to participation in the study. All components of this study were performed in compliance with the Declaration of Helsinki and all procedures were approved by Partners Human Research Committee.
Our study was divided into Part 1 and Part 2. In Part 1, we aimed to correlate the RGE metrics including bER, ΔPO 2 , ΔPCO 2 and ToB with cerebral hemodynamic responses including CBFv and BOLD signal changes to breath hold challenge. We also examined the temporal features and frequency characteristics of RGE metrics and their coherence with cerebral hemodynamic responses. Sixteen out of 17 subjects performed breath hold tasks in the MRI sessions, while 12 of them participated in TCD sessions. In Part 2, we aimed to assess the usefulness of bER in the regional breath hold CVR quantification by comparing the regional CVR maps obtained by regressing BOLD signal changes on bER, P ET CO 2 and ToB under breath hold challenge with regional CVR maps obtained by regressing BOLD signals changes on P ET CO 2 under exogenous CO 2 challenge. Ten out of 17 subjects had additional exogenous CO 2 challenge for comparison. Before we correlated the changes of RGE metrics with CBFv and BOLD signal changes under breath hold challenge, we examined the correlations among the RGE metrics (bER, ΔPCO 2 , ΔPO 2 , P ET CO 2 , P ET O 2 and ToB) acquired in both TCD and MRI sessions.

Part 1: Breath hold challenge
Transcranial Doppler scanning. Before the study of blood flow velocity in intracranial arteries, subject was allowed to rest at least 20-30 minutes for hemodynamic stabilization. The blood pressure measured in the subject was within the normal range [36]. With the subject in a sitting position, a dual probe setting with 2MHz transducers in conjunction with TCD system (Delicate EMS-9U, Shenzhen, China) was used for simultaneous recording of CBFv in the middle cerebral arteries (MCA) on both left and right sides while the subject was performing the breath hold task. Two transducers were attached onto the left and right temporal bone window by velcro. The depth of the Doppler samples was confined to the M1 segment, which is at the main stem of the MCA, for all the subjects.
The timing of the breath hold task was presented visually to the subject by a computer using the software Eprime Professional 2.0 (Psychology Software Tools, Inc., Pittsburgh, USA). A rehearsal session was given to each subject to practice breath hold task. Each subject was instructed via visual cues to perform 6 epochs of 30-second breath hold interleaved with 60-90 seconds of free breathing (Fig 1). They were instructed by visual cues to only hold their breath for as long as they could during the 30-second period. Multiple epochs of breath holding followed by free breathing increased the samples for quantitative analysis. The total duration of the breath hold protocol lasted 10 minutes.
Physiological changes including PCO 2 , PO 2 , electrocardiogram (ECG) and peripheral blood pressure were measured simultaneously with TCD acquisition. A small nasal tubing was placed at the subject's nostril to sample PCO 2 and PO 2 via gas analyzers (Capstar-100, Oxystar-100, CWE, Inc., PA, USA) after calibrating to the barometric pressure of the day of TCD scanning and correcting for vapor pressure. Peripheral blood pressure was continuously measured with Finometer MIDI (Finapres Medical Systems B.V., Netherlands). All the TCD and physiological measurements were synchronized using trigger signals from E-prime. CBFv time series and physiological recordings were stored for offline data analysis.
MRI acquisition. MRI brain scanning was performed on a 3-Tesla scanner (Siemens Medical, Erlangen, Germany). The head was immobilized in a standard head coil with foam pads. The following whole brain MRI datasets were acquired on each subject: 1) standard high-resolution sagittal images acquired with volumetric T1-weighted 3D-MEMPRAGE sequence (TR = 2530ms, TE = 1.74ms/3.6ms/5.46ms/7.32ms, flip angle = 7˚, FOV = 256×256mm, matrix = 256×256, slice thickness = 1mm); 2) BOLD-fMRI images acquired with gradient-echo echo planar imaging (EPI) sequence (TR = 1450ms, TE = 30ms, flip angle = 90˚, FOV = 220×220mm, matrix = 64×64, slice thickness = 5mm, slice gap = 1mm) while the subject was performing the breath hold task. The breath hold task and the physiological set-up used for gas sampling in MRI session were the same as those used in TCD sessions. The gas analyzers were again calibrated to the barometric pressure of the day of MRI scanning and corrected for vapor pressure. ECG was measured using Siemens physiological monitoring unit. Physiological changes including PCO 2 , PO 2 and ECG were measured simultaneously with MRI acquisition. All the physiological measurements were synchronized using trigger signals from the MRI scanner. BOLD-fMRI images and physiological recordings were stored for offline data analysis.

Part 2: Exogenous CO 2 challenge
Ten out of 16 subjects had additional exogenous CO 2 challenge in the MRI session. Given that there is significant inter-individual variance in resting P ET CO 2 [37], resting P ET CO 2 was assessed in those subject via calibrated capnograph before the exogenous CO 2 challenge. Subject wore nose-clip and breathed through a mouth-piece on an MRI-compatible circuit designed to maintain the P ET CO 2 within ± 1-2 mmHg of target P ET CO 2 [38,39]. The fraction of inspired carbon dioxide was adjusted to produce steady-state conditions of normocapnia and mild hypercapnia (4-8 mmHg above the subject's resting P ET CO 2 ). The CO 2 challenge paradigm consisted of 2 consecutive phases (normocapnia and mild hypercapnia) repeating 6 times with 3 epochs of 4 mmHg increase and 3 epochs of 8 mmHg increase of P ET CO 2 (Fig 1). The normocapnia phase lasted 60-90 seconds, while the mild hypercapnia phase lasted 30 seconds. The total duration of the exogenous CO 2 hypercapnic challenge lasted 10 minutes.
When the subject had exogenous CO 2 challenge in MRI session, BOLD-fMRI images were acquired using the same EPI sequence for breath hold challenge. The PCO 2 and PO 2 were sampled through the air filter connected with the mouthpiece and the sampled gases were measured by calibrated gas analyzers. The respiratory flow was measured with respiratory flow head (MTL300L, AdInstruments, Inc., CO, USA) on the breathing circuit via calibrated spirometer (FE141, AdInstruments, Inc., CO, USA). The physiological measurements were synchronized with MRI images using trigger signals from MRI scanner. All the BOLD-fMRI images and physiological recordings were stored for offline data analysis.

Data analysis
Processing of physiological data. The physiological data from both TCD and MRI sessions were analyzed using Matlab R2014a (Mathworks, Inc., Natick, MA, USA). Technical delay of PCO 2 and PO 2 was corrected by cross-correlating the time series of PCO 2 and PO 2 with respiratory phases determined from the artifactual displacement due to chest excursion on ECG time series for the breath hold runs in the TCD sessions, with the respiratory phases from the respiratory bellow for the breath hold runs in the MRI sessions, and with the respiratory flow for the exogenous CO 2 runs.
End inspiration (I) and end expiration (E) were defined on the time series of PO 2 and PCO 2 (S1 Fig). They were verified by the inspiratory and expiratory phases on the respiration time series. The breath-by-breath end-tidal CO 2 (P ET CO 2 ) and end-tidal O 2 (P ET O 2 ) were extracted at the end expiration of PCO 2 and PO 2 time series respectively. Changes of the gas parameters during breath hold periods were interpolated by the values measured immediately before and after the breath hold periods. The duration of breathing cycle, which is known as time of breath (ToB), was derived by subtracting the timing in seconds of the 2 consecutive end expiration markers. Breath-by-breath O 2 -CO 2 exchange ratio (bER) is defined as the ratio of the change in PO 2 (ΔPO 2 = inspired PO 2 -expired PO 2 ) to the change in PCO 2 (ΔPCO 2 = expired PCO 2 -inspired PCO 2 ) measured between end inspiration and end expiration in each respiratory cycle.
Simple correlation analyses were applied on the time series of RGE metrics (bER, ΔPCO 2 , ΔPO 2 , P ET CO 2 , P ET O 2 and ToB) in pairs, which were acquired in both TCD and MRI sessions. Time series of the RGE metrics including bER, ΔPO 2 , ΔPCO 2 and ToB were also used in the CBFv and BOLD data analyses.
Preprocessing of CBFv data. The CBFv data were analyzed using Matlab R2014a (Mathworks, Inc., Natick, MA, USA). A median filter was applied to the data to reduce artifactual spikes. Beat-by-beat systolic peaks and end-diastolic troughs were determined using custom Matlab function and corrected on the graphical user interface incorporated in the function. Systolic peaks and end-diastolic troughs of cardiac cycles on the CBFv time series showing persistent artifacts were excluded in the following analysis. TCD data in both left and right MCAs were acquired on 12 subjects. One of the 12 TCD datasets had persistent artifacts in over onethird of the time series acquired in the LMCA. The CBFv time series in the LMCA of that particular TCD dataset was excluded in further analysis. Time series of mean CBFv was derived by averaging the CBFv over each cardiac cycle. In order to reduce the large inter-individual variations of absolute blood flow velocities [40,41] and to remove the dependence of insonation angle [42], the percent change of CBFv (ΔCBFv) of the left and right MCAs relative to baseline value was derived. The mean CBFv for a period of 30 seconds at the beginning of the time series was chosen as the baseline because the subject was in resting state and the CBFv acquired within this period had the least effect from the respiratory challenges.
Preprocessing of BOLD-fMRI data. All the BOLD-fMRI data were imported into the software Analysis of Functional NeuroImage (AFNI) [43] (National Institute of Mental Health, http://afni.nimh.nih.gov) for time-shift correction, motion correction, normalization and detrending. The first 12 volumes in the first 12 time points of each functional dataset, collected before equilibrium magnetization was reached, were discarded. Each functional dataset was corrected for slice timing, motion-corrected and co-registered to the first image of the first functional dataset using three-dimensional volume registration. It was then normalized to its mean intensity value across the time-series. Voxels located within the ventricles and outside the brain defined in the parcellated brain volume using FreeSurfer [44] (MGH/MIT/HMS Athinoula A. Martinos Center for Biomedial Imaging, Boston, http://surfer.nmr.mgh.harvard. edu) were excluded from the following analyses of functional images. The time-series of each voxel in the normalized functional dataset was detrended with the 5 th order of polynomials to remove the low drift frequency. Individual brain volumes with time series of percent BOLD signal changes (ΔBOLD) were derived.

Part 1
Correlation analysis between cerebral hemodynamic responses and RGE metrics. The cerebral hemodynamic responses (ΔCBFv and ΔBOLD) were separately correlated with the RGE metrics including bER, ΔPCO 2 , ΔPO 2 and ToB. The correlation indicated by Pearson's correlation coefficient was considered significant at p<0.05. Fisher Z-transformation was used to transform Pearson's correlation coefficients to Fisher Z scores for group analysis. Paired ttests were used to compare the Fisher Z scores representing the correlation between cerebral hemodynamic responses and bER with those indicating the correlation between cerebral hemodynamic responses and RGE metrics other than bER. Differences were considered to be significant at p<0.05.
Dynamic analysis of coherence between cerebral hemodynamic responses and RGE metrics as function of time and frequency. Wavelet transform coherence is a method for analyzing the coherence and phase lag between two time series as a function of both time and frequency [45][46][47]. It is therefore well suited to investigating non-stationary changes in coupling between the time series of cerebral hemodynamic responses (ΔCBFv and ΔBOLD) and the time series of RGE metrics including bER, ΔPCO 2 , ΔPO 2 and ToB, as well as the phase lag of the cerebral hemodynamic responses to each of the metrics. The time series of ΔBOLD were extracted from left gray matter (LGM), right gray matter (RGM), left white matter (LWM) and right white matter (RWM). We used the Matlab wavelet cross-spectrum toolbox developed by Grinsted et al. [46]. Squared wavelet coherence between time series of RGE metrics and cerebral hemodynamic responses (ΔCBFv and ΔBOLD) was separately plotted with x-axis as time and y-axis as scale which had been converted to its equivalent Fourier period. An example of squared wavelet coherence between bER and ΔCBFv in right MCA from a representative subject under breath hold challenge is shown in S2 Fig. The magnitude of wavelet transform coherence ranged between 0 and 1 that can be conceptualized as a localized correlation coefficient in time and frequency space [46]. The phase angle between the two time series at particular samples of the time-frequency plane is indicated by an arrow: a rightward pointing arrow indicates that the time series are in phase, or positively correlation (ϕ = 0); a leftward pointing arrow indicates anticorrelation (ϕ = π), and the downward and upward pointing arrows indicate phase angles of π/2 and -π/2 relative to ϕ = 0, respectively. Areas inside the 'cone of influence', which are locations in the time-frequency plane where edge effects give rise to lower confidence in the computed values, are shown in faded color outside of the conical contour. The statistical significance level of the wavelet coherence is estimated using Monte Carlo method and the 5% significance level against red noise is shown as thick contour in the squared wavelet coherence plot. The wavelet coherence magnitudes and phases bounded by thick contour outside the cone of influence are considered significant.
Time-averaged coherence is defined as the total significant coherence at each scale of Fourier periods (converted into frequency) where the wavelet coherence magnitude exceeds 95% significance level, normalized by the maximum possible coherence outside the cone of influence at that particular scale (S2 Fig). It is interpreted in the similar way as the coherence in the transfer function analysis which has been used in cerebral autoregulation study [48]. The mean time-averaged coherence at the phase lags of 0±π/2 and π±π/2 were plotted for all subjects included in TCD or MRI sessions to explore the Fourier periods/frequency bandwidths that oscillations of cerebral hemodynamic responses (ΔCBFv and ΔBOLD) were in synchrony with each physiological time series of bER, ΔPO 2 , ΔPCO 2 and ToB when the subjects were performing breath hold task.

Part 2
Linear regression for CVR quantification. Linear regression analysis was used to derive CVR from the time series of ΔBOLD and vasoactive stimulus when the subject was under breath hold and exogenous CO 2 challenges. The time series of each of the commonly used vasoactive stimuli (bER, P ET CO 2 and ToB) served as a regressor in a separate linear regression analysis. CVR was defined as the percent BOLD signal changes per unit change of the vasoactive stimulus. Therefore CVR was quantified by the coefficient of regression, i.e. the slope.
The statistical parametric maps for individual subjects were cluster-corrected using a threshold estimated with Monte Carlo simulation algorithm. Individual subject brain volume with CVR magnitude was registered onto each subject's anatomical scan and transformed to the standardized space of Talairach and Tournoux [49]. In order to protect against type I error, individual voxel probability threshold of p<0.005 was held to correct the overall significance level to α<0.05. Monte Carlo simulation was used to correct for multiple comparisons [50]. Based upon a Monte Carlo simulation with 2000 iteration processed with ClustSim program [51], it was estimated that a 476mm 3 contiguous volume would provide the significance level α<0.05, which met the overall corrected threshold of p<0.05.
Group comparison of CVR between breath hold and exogenous CO 2 challenges. For each subject who participated in both breath hold and exogenous CO 2 MRI scanning, CVR values were derived from regressing ΔBOLD on bER (CVR BH-bER ), ΔBOLD on P ET CO 2 (CVR BH-PETCO2 ) and ΔBOLD on ToB (CVR BH-ToB ) when the subjects were under breath hold challenge. CVR BH-bER , CVR BH-PETCO2 and CVR BH-ToB values were separately averaged in each of the 160 brain regions parcellated by the software FreeSurfer. Using the same method, CVR values during exogenous CO 2 challenge were obtained by regressing ΔBOLD on P ET CO 2 (CVR CO2-PETCO2 ). CVR CO2-PETCO2 values were averaged in each of the 160 brain regions. To study the CVR changes in group, one-sample t-tests were applied onto the brain volumes with regional CVR BH-bER , CVR BH-PETCO2 , CVR BH-ToB and CVR CO2-PETCO2 . Differences were considered significant at false discovery rate adjusted p fdr <0.05.
The physiological mechanisms underlying breath holding and exogenous CO 2 hypercapnia are potentially different. To study the precision of vasoactive stimulus in predicting the regional cerebral hemodynamic responses, percentage of voxels in each region having significant CVR changes that survived at cluster-corrected p<0.05 in individual subject analysis, in short vCVR, were calculated. Individual subject brain volumes with regional vCVR due to P ET CO 2 , (vCVR BH-PETCO2 ), ToB (vCVR BH-ToB ) and bER (vCVR BH-bER ) were obtained in breath hold MRI scanning, while those with regional vCVR due to P ET CO 2 (vCVR CO2-PETCO2 ) were obtained in exogenous CO 2 MRI scanning. One sample t-tests were again applied for group analysis. Differences were considered significant at p fdr <0.05.
To study the usefulness of P ET CO 2 , ToB and bER in CVR quantification under breath hold challenge, using vCVR CO2-PETCO2 as reference, paired t-tests were applied to compare the regional brain maps of vCVR CO2-PETCO2 under exogenous CO 2 challenge with those of vCVR BH-PETCO2 , vCVR BH-ToB and vCVR BH-bER in groups. Differences were considered significant at p fdr <0.05.

Part 1
Subject demographics were shown in Table 1. The P ET CO 2 , P ET O 2 , ΔPCO 2 , ΔPO 2 and bER measured in the baseline periods, i.e. the period before the first breath hold epoch, in the TCD and MRI sessions were summarized in Tables 2 and 3 respectively. The breath hold duration, the changes of ΔPCO 2 , ΔPO 2 , bER, ΔCBFv and ΔBOLD averaged over the 6 breath hold epochs from the onset to the end of each breath hold epoch were also included. Most of the subjects were able to hold their breaths for 30 seconds on average. The averaged changes of ΔPO 2 from the onset to the end of the breath hold epochs were almost 3 to 4 folds larger than those of ΔPCO 2 .
Correlation among RGE metrics. The correlations among the RGE metrics (bER, ΔPCO 2 , ΔPO 2 , P ET CO 2 , P ET O 2 , and ToB) in TCD and MRI sessions are shown in S1 Fig. Strong and positive correlation was found between bER and ΔPO 2 consistently in all subjects (Pearson's r, 0.70-0.98, p<0.001), while the correlation between bER and ΔPCO 2 varied from weak to moderate (Pearson's r, 0.07-0.75) (S1 Table). Moderate to strong correlation was observed between ΔPO 2 and ΔPCO 2 (Pearson's r, 0.44-0.88), and between P ET O 2 and P ET CO 2 (Pearson's r, 0.46-0.91). Such correlations suggest that ΔPO 2 and ΔPCO 2 (or P ET O 2 and P ET CO 2 ) are not necessarily redundant. Similar observation of non-redundancy between P ET O 2 and P ET CO 2 (correlation coefficients, 0.25-0.93) during spontaneous breathing was also reported by Lenfant et al. [30]. The different ranges of correlation strength observed between ΔPO 2 and ΔPCO 2 in the TCD and MRI sessions may be due to the interaction of the subjects with the environment. Subjects were sitting in an open and quiet environment in the TCD sessions while they were in supine position in a noisy MRI scanner bore. Many previous studies reported that a change from supine posture to sitting upright was associated with a redistribution of both blood flow and ventilation in the lungs, which affected the arterial PO 2 [52][53][54][55]. While it is interesting to observe a difference in the interaction between subjects and environment, the details of such mechanisms are outside the scope of our current study.
Correlation between cerebral hemodynamic responses and RGE metrics. Among the four RGE metrics of bER, ΔPO 2 , ΔPCO 2 and ToB, bER is the only parameter that consistently showed strong correlation with the ΔCBFv measured in left and right MCAs (Pearson's r, 0.40-0.89, p<0.001) (Fig 3A). The bER also correlated with ΔBOLD extracted from gray (LGM and RGM) and white matter (LWM and RWM) under breath hold challenge (Pearson's r, 0.21-0.83, p<0.05), except one subject in the left brain (Pearson's r, -0.011-0.18, p>0.05) (Fig 3A). The results of correlation analyses between cerebral hemodynamic responses and RGE metrics were summarized in S2 Table. Group comparisons of the correlation between cerebral hemodynamic responses (ΔCBFv and ΔBOLD) and RGE metrics (bER, ΔPO 2 , ΔPCO 2 and ToB) also showed that the correlation between cerebral hemodynamic responses and bER was significantly stronger than those between cerebral hemodynamic responses and the RGE metrics other than bER (S2 Table). Although ToB was not as good as bER in correlating with ΔCBFv or ΔBOLD, it was better than ΔPCO 2 to serve as an indicator of the breath hold periods. Dynamic analysis of coherence between cerebral hemodynamic responses and RGE metrics as function of time and frequency. The mean time-averaged coherence between time series of RGE metrics (bER, ΔPO 2 , ΔPCO 2 and ToB) and cerebral hemodynamic responses (ΔCBFv and ΔBOLD) was found to be significantly stronger between 0.008Hz (1/128 seconds) and 0.03Hz (1/32 seconds) at phase lag of 0±π/2 (S3 Fig). We therefore focused on the distribution of time-averaged coherence between RGE metrics and cerebral hemodynamic responses at the phase lag of 0±π/2 (Fig 3). Among the 4 RGE metrics, the mean time-averaged coherence between bER and cerebral hemodynamic responses at phase lag of 0±π/2 was the strongest at the frequency bandwidths of 0.008-0.03 Hz while that between ΔPCO 2 and cerebral hemodynamic responses was the weakest at the same frequency bandwidths. The mean time-averaged coherence between bER, ΔPO 2 , ToB and cerebral hemodynamic responses at phase lag of 0±π/2 reached 0.6 or above. at the frequency bandwidths of 0.008-0.03 Hz. The correlation and dynamic coherence results suggest that CBFv and BOLD signals oscillated with bER at a broad frequency range of low frequencies when the subjects were performing breath hold task. The differences in the correlation and coherence findings between ΔPO 2 and ΔPCO 2 further suggest that changes of PO 2 and PCO 2 are not simply the inverse of each other.

Part 2
Group comparison of CVR between breath hold and exogenous CO 2 challenges. Fig 4 shows the regional CVR brain maps averaged across 16 subjects who performed breath hold task in MRI sessions. Under breath hold challenge, most of the brain regions showed significant increase in CVR BH-bER and CVR BH-ToB especially in thalamus, insula and putamen, while no significant changes of CVR BH-PETCO2 were observed in the same subject group (Fig 4).
In the comparison of CVR brain maps under breath hold and exogenous CO 2 challenge, a subset of 10 subjects who participated in both challenges were included. Under exogenous CO 2 challenge, most of the brain regions showed increased CVR CO2-PETCO2 in the subject group especially thalamus, insula and putamen (Fig 5A), and the vCVR CO2-PETCO2 which had significant CVR CO2-PETCO2 changes exceeded 80% in most of the brain regions (Fig 5B). For  the same group of subjects under breath hold challenge, increased CVR BH-bER and CVR BH-ToB were found in most of the brain regions, while no significant changes of CVR BH-PETCO2 were shown in most of the brain regions ( Fig 5A). CVR brain maps under breath hold challenge shown in Fig 5A (n = 10) were consistent with those shown in Fig 4 (n = 16). Comparing with the vCVR BH-ToB and vCVR BH-PETCO2 , vCVR BH-bER showed the largest percentage of voxels with significant CVR changes in different brain regions, implying a significantly high precision of bER predicting regional ΔBOLD in breath hold challenge (Fig 5B). The paired comparison between vCVR CO2-PETCO2 and vCVR BH-bER , as well as that between vCVR CO2-PETCO2 and vCVR BH-ToB , did not show significant difference in most of the brain regions, while significant differences were found in most of the brain regions between vCVR CO2-PETCO2 and vCVR BH-PETCO2 (Fig 5C). Such findings suggest that bER in breath hold challenge is more appropriate to be used as vasoactive stimulus than P ET CO 2 in assessing regional CVR under breath hold challenge.

Discussion
Our findings show a strong positive correlation between the cerebral hemodynamic responses and our new breath-by-breath O 2 -CO 2 exchange ratio, in short bER, under brief breath hold challenge. We are the first to show that the dynamic changes in bER robustly characterized CBFv and BOLD responses much better than changes in P ET CO 2 or ToB under breath hold challenge in the very low frequency range of 0.008-0.03Hz. The difference between bER and ΔPCO 2 in coherence with cerebral hemodynamic responses within the frequency range of 0.008-0.03 Hz cannot be attributed to the long periods of the breath hold protocol alone since , and the bottom panel shows the coherence with ΔBOLD in LWM and RWM (n = 16). The mean time-averaged coherence in the frequency bandwidths from 0.008 to 0.25Hz were plotted (thick color lines). Color shaded areas represent standard error of the mean. The mean time-averaged coherence between bER and cerebral hemodynamic responses reached 0.6 or above at the frequency range from 0.008Hz (1/128 seconds) to 0.03Hz (1/32 seconds), while the mean time-averaged coherence between ΔPCO 2 and cerebral hemodynamic responses stayed below 0.4. Comparing with ΔPCO 2 , the mean time-averaged coherence of ΔPO 2 with cerebral hemodynamic responses reached 0.5 or above in the frequency range of 0.008-0.03Hz, which was better than that of ΔPCO 2 . https://doi.org/10.1371/journal.pone.0225915.g003 the influence of the periods of the breath hold protocol should contribute equally to all RGE metrics. During breath holding, we presented the combined effect of both hypoxia and hypercapnia on the cerebral hemodynamic responses measured using TCD and BOLD-fMRI. Given that the concurrent changes of P ET O 2 and P ET CO 2 are in opposite direction and the magnitudes depend on the respiratory phase and volume, ΔPO 2 and ΔPCO 2 were in phase and selected here to more conveniently characterize the breath-by-breath changes during breath holding and spontaneous breathing epochs. bER was selected with the O 2 term as the numerator of the O 2 / CO 2 ratio due to bER's positive temporal relationship with CBF changes. In characterizing regional BOLD signal changes to brief breath hold challenge, bER which took into account the interaction of ΔPO 2 and ΔPCO 2 yielded much better results than what P ET CO 2 and ToB could do as we showed that the brain regions outlined by bER during brief breath hold challenge were comparable with those outlined by the P ET CO 2 during exogenous CO 2 challenge. Group CVR maps (CVR BH-PETCO2 , CVR BH-ToB and CVR BH-bER ) generated by regressing ΔBOLD separately on P ET CO 2 , ToB and bER under breath hold challenges (n = 16). All the CVR maps had been corrected for p fdr <0.05. Maps of CVR BH-ToB and CVR BH-bER were comparable, while P ET CO 2 was not able to characterize the regional ΔBOLD. https://doi.org/10.1371/journal.pone.0225915.g004

Mild hypoxia enhances sensitivity of CBF changes to CO 2
During breath holding periods, the CO 2 release are dependent on O 2 uptake in the closed circuit of systemic circulation created by holding one's breath. The interaction between ΔPO 2 and ΔPCO 2 during breath hold is mainly resulting from the systemic metabolic process and different from the effect of exogenous gas administration which is primarily indicated by an increase in ΔPCO 2 . In the current study, ΔPO 2 could go from 11 to 51 mmHg and ΔPCO 2 could go from 2 to 10 mmHg at the end of 30 seconds of breath holding in the TCD and MRI sessions (Tables 2 and 3). This relatively modest level of change in ΔPO 2 had been reported to be able to induce a progressive increase of CBF in the presence of mild hypercapnia [27,28]. Our findings of ΔPO 2 vs. ΔPCO 2 support the synergistic effect of hypoxia and hypercapnia on CBF change under breath hold challenge which are in parallel with the increased chemoreceptor activities in the presence of both mild hypoxia and mild hypercapnia reported by several research teams on animal models [15,17,18,56]. Sensitivity of peripheral chemoreceptor at carotid body to arterial PCO 2 was found to increase when there was a decrease in the level of arterial PO 2 in these studies [15,17,18,56]. The chemoreceptor activity in response to arterial PO 2 took place at all levels of PO 2 ranging from hypoxia to normoxia and then up to 190mmHg or even higher in the range of hyperoxia. The chemoreceptor response curve to PO 2 was similar to a hyperbola with chemoreceptor activity rising faster under hypoxia (below 95mmHg) than under normoxia or hyperoxia [56]. Modulation of peripheral chemoreceptor activities may be expected leading to modulation of CBF while their physiological mechanisms remain to be clarified. . The percentage of voxels with significant CVR changes under breath hold challenge were found in an increasing order of vCVR BH-PETCO2 , vCVR BH-ToB and vCVR BH-bER . CVR BH-bER and vCVR BH-bER showed large resemblance with CVR CO2-PETCO2 and vCVR CO2-PETCO2 respectively. (C) Paired comparisons of vCVR maps showed that vCVR CO2-PETCO2 maps under exogenous CO 2 challenge were not significantly different from vCVR BH-bER maps under breath hold challenge. Group map of vCVR BH-PETCO2 showed that ΔBOLD in only small number of voxels within the brain regions associated with the P ET CO 2 changes under breath hold challenge, resulting in significant differences found in multiple brain regions in the paired comparison between vCVR CO2-PETCO2 and vCVR BH-PETCO2 . https://doi.org/10.1371/journal.pone.0225915.g005

Wavelet transform coherence (WTC) analysis showed strong coherence between bER and cerebral hemodynamic responses (ΔCBFv and ΔBOLD) under breath hold challenge
We used wavelet transform coherence analysis to examine the temporal features and frequency characteristics of bER, ΔPO 2 , ΔPCO 2 and ToB, and their coherence with CBFv and BOLD signal changes under breath hold challenge. Using WTC, we showed that the coherence between cerebral hemodynamic responses (both ΔCBFv and ΔBOLD) and change in bER was much stronger than that between cerebral hemodynamic responses and change in ΔPCO 2 in a wide frequency range of 0.008-0.03Hz (Fig 3). The meaning of the range of 0.008-0.03Hz is best appreciated by examining the WTC findings which showed that bER and ΔPCO 2 were very different in the frequency distribution of their coherence with cerebral hemodynamic responses. The special characteristics of the frequency range of 0.008-0.03Hz need to be considered in the context of each individual RGE metric (ΔPO 2 , ΔPCO 2 and bER) separately.
Comparing with ΔPCO 2 , the stronger coherence found between bER and cerebral hemodynamic responses in low frequency range of 0.008-0.03Hz may be associated with low frequency physiological processes in the brain that are grouped in B-wave frequency bandwidth. B-waves with a period of 0.5 to 2 minutes have been reported to be related to autoregulation of microvasculature and spontaneous rhythmic oscillations in intracranial pressure [57,58]. Changes in systemic PO 2 and PCO 2 may trigger the underlying mechanisms which involve the interaction of central and peripheral respiratory chemoreceptors as well as autonomic system to regulate CBF. B-waves may be associated in such a process through the contractile properties of pericytes or vascular smooth muscle cells to alter vascular diameter and ensure the delivery of O 2 and other circulating metabolites [59][60][61][62]. As autoregulation of CBF and oxygen delivery are the outcomes in this part of vital homeostatic process [17], it is reasonable that the stronger coupling was found between bER and cerebral hemodynamic responses, and between ΔPO 2 and cerebral hemodynamic responses.
In summary, among the four RGE metrics of bER, ΔPO 2 , ΔPCO 2 and ToB, the cerebral hemodynamic responses showed the strongest correlation and dynamic coherence with bER, followed by ΔPO 2 , ToB, and ΔPCO 2 . A wide range of correlations found between ΔPO 2 and ΔPCO 2 (and between P ET O 2 and P ET CO 2 ) in the subject group, as well as the differences in the correlation in the time domain and coherence findings in wavelet transform coherence analysis between ΔPO 2 and ΔPCO 2 indicate that PO 2 and PCO 2 are not simply inverse of each other. The strong correlation between bER and ΔPO 2 (S1 Fig) indicates that bER is predominantly affected by ΔPO 2 . ToB yielded superior correlation result than ΔPCO 2 because ToB indirectly takes into account the duration for both hypoxemia and hypercapnia which both elevate cerebral hemodynamic responses without being affected by the depth of breathing. bER is the most accurate in characterizing cerebral hemodynamic responses under breath hold challenge because it directly takes into account the effect of both ΔPO 2 and ΔPCO 2 .

CVR map from breath hold challenge resembled that from exogenous CO 2 challenge
Our CVR findings showed that there was a large resemblance between regional CVR BH-bER and CVR CO2-PETCO2 , even though the underlying physiological mechanisms for the CBF increase are potentially different (Fig 5). From spontaneous breathing to breath holding periods, the bER, ΔPO 2 and ΔPCO 2 measured were due to endogenous changes of respiratory gases involved in the process of systemic metabolism. Previous studies reported that apneainduced hypoxia and hypercapnia caused chemoreceptor-mediated central vasodilation and concurrent peripheral vasoconstriction to conserve oxygen delivery to the brain [25], leading to an increase in CBF and a decrease in peripheral oxygen saturation [26]. Holding breath is different from exogenous CO 2 administration where ΔPO 2 and ΔPCO 2 from normocapnic to hypercapnic epochs depended on the externally administered gas mixture and the increase in cerebral hemodynamic responses was mainly due to hypercapnia. One explanation for the large resemblance of the CVR maps between the breath hold and exogenous CO 2 challenges is that bER as an optimal regressor under breath hold challenge is able to characterize the increase in cerebral hemodynamic responses in most of the brain regions. Comparing the brain maps of vCVR BH-bER with those of vCVR BH-ToB and vCVR BH-PETCO2 , more voxels showed significant association between ΔBOLD and bER. For CVR maps, CVR BH-ToB resembled CVR CO2-PETCO2 more than CVR BH-PETCO2 did because ToB is a good time indicator of breath hold epochs as shown in Fig 2. Breath hold challenge, as shown in our findings in Part 2, was able to offer significant regional CVR quantification, as long as the proper regressor was used. Since bER is a measure of the dynamic change of RGE which is related to ratio of the change of partial pressure in O 2 uptake to that in CO 2 release in the body, mapping CVR to breath hold may offer a novel probe to study the interaction between brain responses and some measures of systemic metabolism. Mild exogenous CO 2 challenge is expected to have little effect on cerebral [63] and presumably systemic metabolism. From a technical approach, breath-holding challenge for CVR assessment is much easier to set up in clinical settings than external CO 2 administration which requires an elaborate gas administration circuit [6]. The breath-holding challenge also allows patients to switch back to normal breathing whenever their physiological limit is reached [64]. In patients who have compromised vasculature with potential risk of acute intracranial hemorrhage, breath-holding may offer an option for CVR study if one has some yet-to-be proven concern over the more powerful physiological stress presented by rapid increase of externally administered CO 2 .
Why was the ratio of ΔPO 2 to ΔPCO 2 used, and not the product?
The success of bER in better characterizing dynamic CBFv and BOLD signal changes under breath hold challenge is closely related to the fact that bER is a ratio which factors out at the same time effects of ventilatory volume fluctuations [32] common to both ΔPO 2 and ΔPCO 2 . Given that both hypoxia and hypercapnia induce the increase in cerebral hemodynamic responses, one may wonder whether the product of ΔPO 2 and ΔPCO 2 would be better than bER for the evaluation of change in cerebral hemodynamic responses. Our answer is no. First, the product of ΔPO 2 and ΔPCO 2 would exacerbate the contribution of fluctuations from ventilatory volume. Secondly, unlike the ratio which has long been used to properly describe RGE since at least 1913 [65], the product does not have an established physiological meaning. Actually, another way to look at bER is that it describes the change of ΔPO 2 per unit change of ΔPCO 2 , so the ratio provides a way to quantify how O 2 and CO 2 can work together to interact with cerebral hemodynamic responses.
In a related framework, the 'stimulus index' (SI; P ET CO 2 /P ET O 2 ) developed by Bruce et al. [66] for breath hold study is different from bER, namely ΔPO 2 /ΔPCO 2 , in its interaction with CBF. Both Bruce et al. (using SI) and our team (using bER) agree upon the influence of O 2 and CO 2 on CBFv change in breath hold. But there are several major differences between SI and bER.
One major property that contributes to the differences between bER and SI is the direction/ phase of oscillations of gas measurements. As shown in Fig 2, the time courses of P ET O 2 and P ET CO 2 oscillate out of phase, while ΔPO 2 oscillates in phase with ΔPCO 2 . With ΔPO 2 being in phase with ΔPCO 2 , shared but less interesting physiological signals (e.g. respiratory fluctuation) would be largely factored out in the ratio ΔPO 2 /ΔPCO 2 (bER). With P ET O 2 and P ET CO 2 being out of phase, those less interesting physiological signals do not get factored out by P ET CO 2 /P ET O 2 (SI) but can in fact be exaggerated in amplitude, as indicated by two-headed arrows in red in Figs 2 and S4. Hence those physiological fluctuations that may be less relevant with the direct interaction between cerebral hemodynamic responses and O 2 or CO 2 would be larger for SI than for bER. We believe that is one of the explanations for why bER is superior to SI in its correlation with CBF ( S4 Fig). Another difference is that Bruce et al. [66] studied the correlation between CBF and the ratio P ET CO 2 /P ET O 2 within a single prolonged period (~93 seconds) of breath holding, i.e. from the onset to the end of breath hold only. In our study, we evaluated the relationship between CBF and bER throughout the duration of six breath holding epochs (~30 seconds) and free breathing periods (~60-90 seconds). We took into account the long delayed cerebrovascular response to breath hold as CBF continues to increase shortly after the breath holding period before it slowly returns to baseline during the free breathing period. Our team reported this delayed CBF response in our other TCD study on breath hold in 2009 [67] where we showed that the time duration of CBFv response to breath hold could last twice as long as the time duration of breath hold [67]. Another team also reported the same delayed CBF response to breath hold in an MRI study [10]. That is one reason why our use of bER is expected to be superior to the use of SI developed by Bruce et al. [66] in characterizing the CBF changes in response to breath hold. The stronger positive correlation of CBFv with bER than with SI suggests that there is synergistic behavior of O 2 and CO 2 in bER that goes beyond what SI does.

Conclusion
Independent of our current knowledge to completely clarify why bER offers the strongest association with cerebral hemodynamic responses to breath holding, we succeeded in showing that bER was superior to ΔPCO 2 or P ET CO 2 to characterize cerebral hemodynamic responses under breath hold challenge for the CVR evaluation. RGE metrics of both ΔPO 2 and ΔPCO 2 should always be acquired for CVR evaluation instead of acquiring P ET CO 2 data alone. In addition to offering alternative approach of CVR evaluation for patients who are not eligible for exogenous CO 2 challenge, the association between bER and cerebral hemodynamic responses also provides a novel insight in the study of brain-body interaction. Future studies would be required to clarify the underlying mechanisms for the relationship between dynamic bER and cerebral hemodynamic response to breath holding. Studies to quantify the relationship between changes in cerebral hemodynamic responses during breath holding and changes in bER in a large cohort of subjects and patients would be helpful to explore the effects on CVR by various disorders including respiratory or cerebral diseases with neurovascular deficits.

S1 Fig. Definition of end inspiration and end expiration on the time series of RGE metrics and the correlations among breath-by breath RGE matrices. (A)
A segment of 80-second time series of ΔCBFv in left MCA and physiological changes including breath-by-breath bER, ΔPO 2 , ΔPCO 2 , P ET O 2 and P ET CO 2 measured by gas analyzers and respiration time series (Resp) measured by respiratory bellow in a representative subject under breath hold challenge in TCD session. Open circles represent end expiration while closed circles represent end inspiration in resting phase or onset of expiration at the end of breath hold epoch. Positive phases with deflection above zero on the respiration time series represent inspiration and negative phases with deflection below zero represent expiration. The inspiratory and expiratory phases of each respiratory cycle on the time series of P ET O 2 and P ET CO 2 are verified by those on respiration time series. The timing for open (end expiration) and closed (end inspiration) circles in green is the same as those in red and blue. (B) Correlations among breath-by breath respiratory matrices (bER, ΔPO 2 , ΔPCO 2 , ToB, P ET O 2 and P ET CO 2 ) in all subjects who participated in TCD sessions (n = 12), and (C) those who participated in MRI sessions (n = 16). Each gray circle represents the Pearson's correlation coefficient from the correlation analysis of the time series of parameter pair shown on x-axis for each subject. The thick middle horizontal line, the box and the vertical rod represent the mean, 95% confidence interval and standard deviation of the group data respectively. The time series of bER had stronger correlation with that of ΔPO 2 than ΔPCO 2 , although both ΔPO 2 and ΔPCO 2 contributed to changes of bER. The correlation coefficients from ΔPO 2 vs ΔPCO 2 varied from 0.6 to 0.9 in TCD sessions and from 0.4 to 0.9 in MRI sessions, suggesting that ΔPO 2 and ΔPCO 2 are not necessarily redundant. The difference in the ranges of correlation strength found between TCD and MRI sessions may be due to the difference in posture of the subjects, where the subjects were in erect seated position in TCD sessions and they were in supine position in MRI sessions. The squared wavelet coherence between these two time series. Squared wavelet coherence is plotted with x-axis as time and y-axis as scale which has been converted to its equivalent Fourier period. The magnitude of wavelet transform coherence ranges between 0 and 1, where warmer color represents stronger coherence and cooler color represents weaker coherence. Areas inside the 'cone of influence', which are locations in the time-frequency plane where edge effects give rise to lower confidence in the computed values, are shown in faded color outside of the conical contour. The statistical significance level of the wavelet coherence is estimated using Monte Carlo methods and the 5% significance level against red noise is shown as thick contour. The phase angle between the two time series at particular samples of the time-frequency plane is indicated by an arrow (rightward pointing arrows indicate that the time series are in phase or positively correlation, leftward pointing arrows indicate anticorrelation and the downward pointing arrows indicate phase angles of π/ 2). There are four different ranges of phase lags: 0+π/2, 0-π/2, π-π/2, and π+π/2. (C) Timeaveraged coherences at four different phase lags of 0+π/2, 0-π/2, π-π/2, and π+π/2. At each phase lag range, time-averaged coherence was defined as the total significant coherence at each scale where the wavelet coherence magnitude exceeded 95% significance level, normalized by the maximum possible coherence outside the cone of influence, i.e. inside the conical contour, at that particular scale and phase lag range. (TIF)

S4 Fig. Time series of cerebral hemodynamic responses, bER and P ET CO 2 /P ET O 2 in a representative subject under breath hold challenge.
Time series of ΔCBFv in left MCA, bER and P ET CO 2 /P ET O 2 in the same representative subject in Fig 2A under breath hold challenge in TCD session. Shaded areas represent breath hold periods. The time series of bER followed closely to the ΔCBFv changes, while P ET CO 2 /P ET O 2 did not follow ΔCBFv changes some time during the challenge as indicated by two-headed arrow in red. In the time period between 400 and 580 seconds (as indicated by two-headed arrow) when the subject had shallow breathing, the amplitude of P ET CO 2 /P ET O 2 decreased significantly in comparison with that of bER. This may be attributed to the different property of gas measurements where the time series of ΔPO 2 and ΔPCO 2 oscillated in phase while those of P ET O 2 and P ET CO 2 oscillated out of phase ( Fig  2A). P ET CO 2 was decreased and P ET O 2 was increased by shallow breathing, resulting in a significant decrease in P ET CO 2 /P ET O 2 . (TIF) S1 Table. Correlation among RGE metrics. Strength of correlation indicated by Pearson's correlation coefficients among respiratory metrics including bER, ΔPO 2 , ΔPCO 2 , ToB, P ET O 2 and P ET CO 2 in all subjects who participated in TCD sessions (n = 12), and those who participated in MRI sessions (n = 16). The time series of bER had stronger correlation with that of ΔPO 2 than ΔPCO 2 , although both ΔPO 2 and ΔPCO 2 contributed to changes of bER. The correlation coefficients from ΔPO 2 vs ΔPCO 2 varied from 0.6 to 0.9 in TCD sessions and from 0.4 to 0.9 in MRI sessions, suggesting that ΔPO 2 than ΔPCO 2 are not necessarily redundant. (DOCX) S2 Table. A. Correlation between RGE metrics and ΔCBFv in TCD sessions. Strength of correlation indicated by Pearson's correlation coefficients between ΔCBFv and RGE metrics including bER, ΔPO2, ΔPCO2 and ToB (n = 12). Numbers in brackets next to Pearson's correlation coefficients indicate p values from individual correlation analyses. The bottom row shows the mean values of Fisher Z scores transformed from Pearson's correlation coefficients in groups. Numbers in brackets next to mean Fisher Z scores indicate p values in the paired comparisons. The correlation between ΔCBFv and bER was significantly larger than those of the correlation between ΔCBFv and the other respiratory metrics in the paired comparisons (p<0.001). bER is the only parameter that consistently showed significantly high correlation with the ΔCBFv measured in LMCA and RMCA. B. Correlation between RGE metrics and ΔBOLD in MRI sessions. Strength of correlation indicated by Pearson's correlation coefficients between ΔBOLD and RGE metrics including bER, ΔPO2, ΔPCO2 and ToB (n = 16). Numbers in brackets next to Pearson's correlation coefficients indicate p values from individual correlation analyses. The bottom row shows the mean values of Fisher Z scores transformed from Pearson's correlation coefficients in groups. Numbers in brackets next to mean Fisher Z scores indicate p values in paired comparisons. The correlation between ΔBOLD and bER was significantly larger than those of the correlation between ΔBOLD and the other respiratory metrics in the paired comparisons (p<0.001). bER is the only parameter that consistently showed significantly high correlation with the ΔBOLD measured in LGM, RGM, LWM and RWM.  Kwong.