Dynamic brain-body coupling of breath-by-breath O2-CO2 exchange ratio with resting state cerebral hemodynamic fluctuations

Background The origin of low frequency cerebral hemodynamic fluctuations (CHF) in the resting state remains unknown. Breath-by breath O2-CO2 exchange ratio (bER) has been reported to correlate with the cerebrovascular response to brief breath hold challenge at the frequency range of 0.008–0.03Hz in healthy adults. bER is defined as the ratio of the change in the partial pressure of oxygen (ΔPO2) to that of carbon dioxide (ΔPCO2) between end inspiration and end expiration. In this study, we aimed to investigate the contribution of respiratory gas exchange (RGE) metrics (bER, ΔPO2 and ΔPCO2) to low frequency CHF during spontaneous breathing. Methods Twenty-two healthy adults were included. We used transcranial Doppler sonography to evaluate CHF by measuring the changes in cerebral blood flow velocity (ΔCBFv) in bilateral middle cerebral arteries. The regional CHF were mapped with blood oxygenation level dependent (ΔBOLD) signal changes using functional magnetic resonance imaging. Temporal features and frequency characteristics of RGE metrics during spontaneous breathing were examined, and the simultaneous measurements of RGE metrics and CHF (ΔCBFv and ΔBOLD) were studied for their correlation. Results We found that the time courses of ΔPO2 and ΔPCO2 were interdependent but not redundant. The oscillations of RGE metrics were coherent with resting state CHF at the frequency range of 0.008–0.03Hz. Both bER and ΔPO2 were superior to ΔPCO2 in association with CHF while CHF could correlate more strongly with bER than with ΔPO2 in some brain regions. Brain regions with the strongest coupling between bER and ΔBOLD overlapped with many areas of default mode network including precuneus and posterior cingulate. Conclusion Although the physiological mechanisms underlying the strong correlation between bER and CHF are unclear, our findings suggest the contribution of bER to low frequency resting state CHF, providing a novel insight of brain-body interaction via CHF and oscillations of RGE metrics.


Introduction
Low frequency components (below 0.05 Hz) in the cerebral hemodynamic fluctuations (CHF) are used to characterize functional connectivity of cerebral resting state networks including the default mode network (DMN) [1,2]. The origin of these low frequencies of CHF is still not clear, and many physiological candidates of cardiac and respiratory origin had been proposed contributing to this frequency bandwidth [3][4][5][6]. They include fluctuations related to respiratory variation (~0.03Hz) [3,4], heart rate variability (0.05-0.15Hz) [5] and end-tidal carbon dioxide fluctuations (0-0.05Hz) [6].
Fluctuations of respiratory gas exchange (RGE) are of particular interest because RGE is a physiological process of removing carbon dioxide (CO 2 ) as a metabolic waste from the blood and replenishing oxygen (O 2 ) for energy consumption. Breath-by-breath RGE metrics include parameters like end-tidal partial pressures of O 2 (P ET O 2 ) and CO 2 (P ET CO 2 ), changes of the partial pressure of O 2 (ΔPO 2 ) and CO 2 (ΔPCO 2 ) [7,8]. According to the concept of homeostasis, the homeostatic regulatory system is formed with arterial PO 2 and PCO 2 as regulated variables, peripheral and central chemoreceptors as sensors, brain stem as the control center, and diaphragm and respiratory muscles as effectors, to optimize systemic blood gases. While the detailed mechanisms between the change in central chemoreceptor activities and the change in cerebral blood flow (CBF) are topics of on-going research [9][10][11], arterial PO 2 and PCO 2 which are sensed by chemoreceptors are maintained within a range (or fluctuate) around the physiological 'set point' (i.e. mean) by the feedback control of diaphragm and respiratory muscles for ventilation in the homeostatic regulatory system. Spontaneous breathing is, therefore, part of a vital homeostatic process to optimize the systemic blood gases which can presumably regulate in turn CBF and O 2 delivery to the brain [8,12]. Such homeostatic fluctuations of arterial PO 2 or PCO 2 may interact synergistically and contribute to fluctuations of CHF during spontaneous breathing at rest. RGE may be related to CHF because oscillations of RGE metrics (P ET O 2 , P ET CO 2 , respiratory exchange ratio) had been reported to be below 0.05Hz in humans [13,14]. In our recent breath hold study, we made an emphasis on the ratio of ΔPO 2 to ΔPCO 2 between end inspiration and end expiration which we named breath-by-breath exchange ratio (bER), to quantify RGE as well [15]. We showed that cerebral hemodynamic responses to brief breath hold challenge were separately coupled with ΔPO 2 , ΔPCO 2 and bER. It is interesting that bER was superior to either ΔPO 2 or ΔPCO 2 alone in coupling with the changes of global/regional cerebral hemodynamic response in our breath hold study. In this manuscript, we hypothesized a brainbody interaction during spontaneous breathing where various RGE metrics would separately contribute to CHF and the coupling of CHF with bER would be the strongest. Our study was different from previous works by Wise et al. [6] and Golestani et al. [16] which assumed that changes of PCO 2 and PO 2 were redundant.
The oscillations of P ET O 2 and P ET CO 2 during spontaneous breathing were shown to be non-redundant by at least three research teams [13,17,18] even though the arterial PO 2 and PCO 2 are expected to be correlated with each other through the natural RGE process. Interdependence but non-redundancy were also shown between P ET O 2 and P ET CO 2 , between ΔPO 2 and ΔPCO 2 , and between PO 2 and PCO 2 in healthy subjects under brief breath hold challenge [15]. Notably, had P ET O 2 /ΔPO 2 /PO 2 and P ET CO 2 /ΔPCO 2 /PCO 2 been redundant during spontaneous breathing, there would be no significant oscillations of respiratory exchange ratio as reported by Lenfant et al. [13]. The non-redundancy between P ET O 2 and P ET CO 2 [13,17,18] is highly significant because it challenges the long-held belief that the oscillations of the partial pressure of O 2 (PO 2 ) and CO 2 (PCO 2 ) are simply mirror image of each other [4], and encourages studying the separate effects of O 2 , CO 2 and their ratio bER on CHF during spontaneous breathing.
In addition to the more commonly used P ET O 2 and P ET CO 2 , we focused more on the contribution of RGE metrics including ΔPO 2 and ΔPCO 2 and bER in this manuscript. We preferred to use the ratio of ΔPO 2 /ΔPCO 2 for the respiratory gas exchange instead of P ET CO 2 / P ET O 2 with two reasons that we have explained with the findings in our recent breath hold study [15]. First, bER is a ratio which factors out the effects of ventilatory volume fluctuations [7] common to both ΔPO 2 and ΔPCO 2 , but not to both P ET CO 2 and P ET O 2 . The time courses of P ET CO 2 and P ET O 2 oscillate out of phase, while ΔPO 2 oscillates in phase with ΔPCO 2 . Second, bER is a breath-by-breath dynamic form related to the steady-state respiratory exchange ratio (RER) described in the alveolar gas equation introduced by Fenn et al. [19]; the timeaveraged value of bER over minutes is therefore mathematically equivalent to the reciprocal of RER. ΔPO 2 and not P ET O 2 is the designated term used in RER from the alveolar air equation [7,19]. ΔPO 2 is used as numerator in bER in this manuscript because of the strong positive correlation found previously between bER and cerebral hemodynamic responses under brief breath hold challenge [15]. RER has been used to evaluate resting systemic metabolic rate [20][21][22]. There are also technical differences between RER used in the literature and bER we used in this study as well as our previous breath hold study [15]. 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. In our study, bER was derived by measuring the inspired and expired gases with a nasal tubing at each breath.
Given the dynamic changes in PO 2 and PCO 2 breath-by-breath in humans, we studied RGE metrics as surrogates of arterial blood gases with the recognition of the difference between respiratory metrics and blood gases. We studied the interaction between RGE metrics and CHF measured with transcranial Doppler sonography (TCD) and functional magnetic resonance imaging (fMRI). TCD has an advantage of acquiring data at high temporal resolution (~100 Hz), eliminating the concern of unwanted high frequency cardiovascular signal being aliased into the low frequency range which has a potential association with the minutelong oscillation of RGE metrics. The cerebral blood flow velocities (CBFv) were measured by TCD in the middle cerebral arteries (MCA) which supply most parts of the brain. The application of TCD allows the subjects to have measurements in an upright seated position, which is the usual position in respiratory physiology studies. Although TCD offers a high temporal resolution to evaluate CHF, it does not provide regional information. Regional mapping of spontaneous CHF was carried out instead with blood oxygen level-dependent (BOLD) signal changes measured by fMRI. We used BOLD-fMRI instead of arterial spin labeling (ASL) in MRI perfusion studies for two major reasons. ASL perfusion-related signal has a low contrast to noise ratio in comparison with the BOLD signal, especially in white matter area. ASL image acquisition at a temporal resolution of 4 seconds also under-samples the spontaneous fluctuations within the respiratory cycle of 4-6 seconds, which hampers the dynamic analyses between CHF and RGE metrics. While CBFv and BOLD signals are not equivalent surrogates for CBF, combining their strengths would be powerful to study the temporal and spatial features in the association between RGE and CHF.
In the present study, we aimed to evaluate the contribution of breath-by-breath RGE metrics to the CHF indicated by CBFv and BOLD signal changes during spontaneous breathing. To address the question of redundancy among RGE metrics during spontaneous breathing, we examined the correlations among bER, ΔPCO 2 , ΔPO 2 , P ET CO 2 and P ET O 2 . In the TCD study, we measured the correlation of CBFv in the MCA with RGE metrics of bER, ΔPO 2 and ΔPCO 2 . We also examined the temporal features and frequency characteristics of these RGE metrics and their coherence with CBFv. In the fMRI study, we mapped the association between regional BOLD signal changes and RGE metrics of bER, ΔPO 2 and ΔPCO 2 . Brain regions showing a significant association between regional BOLD signal changes and RGE metrics were then compared with regions within the default mode network (DMN) outlined by the resting state connectivity analysis with the seed at the left precuneus. DMN encompasses brain regions responsible for the constant background activities at rest, showing higher metabolic and hemodynamic level than the other parts of the cortex [23]. The temporal features and frequency characteristics of these RGE metrics and their coherence with BOLD signal changes in the regions within DMN were also examined. The respiratory metric of respiration volume per unit time (RVT) was also derived to examine if the potential coherence between CHF and RGE metrics were resulting from respiratory variability. The investigation of the association between RGE metrics and CHF would provide new directions to the study of brain-body interaction, in addition to offering a physiological model to characterize the contribution of gas exchange elements in the low frequency resting state fluctuations, With the reference of cerebrovascular reactivity (CVR) to exogenous CO 2 challenge, the potential of using the interaction between CHF and RGE metrics to evaluate CVR to the elements in spontaneous breathing was also discussed.

Participants
Twenty-two volunteers aged from 19 to 48 years (mean age = 30.5 years, SD = 9.1 years, 14 males and 8 females) were included. Eleven of them participated in both TCD and MRI sessions, while the remaining participated in either one of the sessions. Subject demographics were shown in Table 1. 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 Health-Care. All the experimental procedures were explained to the subjects, and signed informed consent was obtained prior to the participation in the study. All components of this study were performed in compliance with the Declaration of Helsinki, and all procedures were approved by the Massachusetts General Hospital IRB of Partners HealthCare.
Our study was divided into two parts: Part I and Part II. In Part I TCD study, we aimed to correlate the RGE metrics including bER, ΔPCO 2 and ΔPO 2 with CBFv in MCAs which serve as CHF in the major cerebral arterial supply during spontaneous breathing. We also examined the temporal features and frequency characteristics of these RGE metrics and their coherence with CBFv. Thirteen subjects participated in TCD sessions. In Part II MRI study, we aimed to map the association between these RGE metrics (bER, ΔPCO 2 and ΔPO 2 ) and BOLD signal changes which serve as regional CHF during spontaneous breathing. Twenty subjects participated in MRI sessions. Ten out of 20 subjects had additional exogenous CO 2 challenge in the MRI sessions. Before we correlated the changes of RGE with CBFv and BOLD signal changes, we examined the correlations among the respiratory metrics (bER, ΔPCO 2 and ΔPO 2 ) acquired in both TCD and MRI sessions.

Part 1: TCD
Transcranial Doppler scanning. Before the study of blood flow velocity in intracranial arteries, the subject was allowed to rest at least 20-30 minutes in an upright seated position for hemodynamic stabilization. The blood pressure measured in the subject was within the normal range of systolic blood pressure below 130 mmHg [24]. With the subject in an upright seated 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 MCA on both left and right sides while the subject was at rest. Two transducers were attached to the left and right temporal bone windows 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.
A white crosshair in the black background was presented visually to the subject by a computer using the software Eprime Professional 2.0 (Psychology Software Tools, Inc., Pittsburgh, USA) to fixate their eye movement. The total duration of the CBFv data acquisition lasted 10 minutes.

PLOS ONE
Dynamic brain-body coupling between resting state hemodynamics and O 2 -CO 2 exchange ratio Physiological changes including PCO 2 , PO 2 , electrocardiogram (ECG) and peripheral blood pressure were measured simultaneously with TCD acquisition. A biological amplifier (FE132, ADInstruments, Inc., CO, USA) and a finometer (Finapres Medical Systems, Netherlands) were used to measure ECG and peripheral blood pressure, respectively. 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 on the day of TCD session and correcting for vapor pressure. Small nasal tubing was chosen to sample PCO 2 and PO 2 because of several reasons. In the current study, we standardized all the physiological set-up in both TCD and MRI sessions. We did not measure PCO 2 and PO 2 from a facemask connected to a closed breathing circuit with non-rebreathing valves because of the space constraint of the MRI head coil, or from the expired gases collected in Douglas bag because of the slow sampling of gases in the units of minutes. Sampling gases with small nasal tubing has an advantage that the subjects breathe normally through the nose. In the gas sampling circuit, the same gas sample volume was used in both CO 2 and O 2 analyzers at the same gas sampling flow rate. The CBFv and physiological measurements were synchronized using trigger signals from E-prime. All the CBFv time series and physiological recordings were stored for offline data analysis.

Part 2: MRI
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 (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 gradientecho echo planar imaging (EPI) sequence (TR = 1450ms, TE = 30ms, flip angle = 90º, FOV = 220×220mm, matrix = 64×64, thickness = 5mm, slice gap = 1mm) while the subject was at rest. The visual presentation of the crosshair, the physiological set-up for the sampling of PCO 2 and PO 2 in the MRI session were the same as those used in the TCD session. The gas analyzers were again calibrated to the barometric pressure on the day of MRI session and corrected for vapor pressure. ECG was measured using a Siemens physiological monitoring unit (Siemens Medical, Erlangen, Germany). Physiological changes including PCO 2 , PO 2 , ECG and respiration were measured simultaneously with MRI acquisition. All the physiological measurements, including those of PCO 2 , PO 2 and ECG were synchronized using trigger signals from the MRI scanner. BOLD-fMRI images and physiological recordings were stored for offline data analysis.
Exogenous CO2 challenge. Ten out of 20 subjects had additional exogenous CO 2 challenge in the MRI sessions. An in-house gas delivery and mixing system comprising of a medical gas mixer in series with a manifold of flow meters was used to mix and deliver gas mixture for exogenous CO 2 challenge. Given that there is significant inter-individual variance in resting P ET CO 2 [25], resting P ET CO 2 was assessed in subjects via calibrated capnograph before the exogenous CO 2 challenge. The 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 [26,27]. 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 . The normocapnic phase lasted 60-90 seconds, while the mild hypercapnic phase lasted 30 seconds. The same paradigm was also used and described in our other study [15]. 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 resting state. 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 the 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 in the TCD sessions, with the respiratory phases from respiratory bellow or respiratory flow measured with respiratory flow head via spirometer in the MRI sessions.
End inspiration (I) and end expiration (E) were defined on the time series of PO 2 and PCO 2 (S1 Fig). The breath-by-breath P ET CO 2 and P ET O 2 were extracted at the end expiration of PCO 2 and PO 2 time series respectively. ΔPO 2 is defined as (inspired PO 2 -expired PO 2 ) and ΔPCO 2 is defined as (expired PCO 2 -inspired PCO 2 ). Breath-by-breath O 2 -CO 2 exchange ratio (bER) is defined as the ratio of ΔPO 2 to ΔPCO 2 measured between end inspiration and end expiration at each breath, whereas RER at steady state from alveolar air equation introduced by Fenn et al. [19] is formulated as the ratio of ΔPCO 2 to ΔPO 2 over minutes. The product of ΔPO 2 and ΔPCO 2 was not used to evaluate the interaction with CHF because the effects of fluctuations due to ventilation would be exacerbated in ΔPO 2 ×ΔPCO 2 (S2 Fig).
Simple correlation analyses were applied to the time series of RGE metrics (bER, ΔPCO 2 , and ΔPO 2 ) in pairs. The correlation was considered significant at p<0.05.

Part 1
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 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 13 subjects. One of the 13 TCD datasets had persistent artifacts in over onethird of the CBFv time series acquired in the LMCA, and another one had persistent artifacts in CBFv data acquired in RMCA. The CBFv time series of those particular TCD runs that showed persistent artifacts were excluded in further analysis, and the CBFv time series acquired in the contralateral MCA that did not show persistent artifacts were retained. Time series of mean CBFv were derived by integrating the CBFv over each cardiac cycle. To reduce the large inter-individual variations of absolute blood flow velocities [28,29] and to remove the dependence of insonation angle [30], the percent change of CBFv (ΔCBFv) relative to baseline value in the left and right MCAs was derived. The mean CBFv for 30 seconds at the beginning of the time series was chosen as the baseline.

Correlation analyses between
ΔCBFv and RGE metrics. The time series of ΔCBFv were correlated with those of bER, ΔPCO 2 and ΔPO 2 separately. The correlation indicated by Pearson's correlation coefficient was considered significant at p<0.05. Fisher's Z-transformation was used to transform Pearson's correlation coefficients to Fisher's z scores for group analysis. Paired t-tests were used to compare the Fisher's z scores representing the correlation between ΔCBFv and bER with those indicating the correlation between ΔCBFv and other physiological parameters besides bER. Differences were considered to be significant at p<0.05. We did not apply a low pass filter to CHF before our correlation analyses, even though applying a low pass filter to CHF is a standard preprocessing step in the resting state analysis [1].
Dynamic analysis of coherence between ΔCBFv and RGE metrics as a function of time and frequency. Wavelet transform coherence (WTC) was employed to demonstrate the dynamic interaction between ΔCBFv and RGE metrics (S3 Fig) in the time-frequency domain. WTC is a method for analyzing the coherence and phase lag between two time series as a function of both time and frequency [31,32]. The temporal and phase information of WTC has been used to map the dynamic connectivity of the brain regions related to heart rate changes [33]. It is therefore well suited to investigate the dynamic changes in the coupling between the time series of ΔCBFv and RGE metrics including bER, ΔPCO 2 , and ΔPO 2 , as well as the phase lag of ΔCBFv to bER, ΔPCO 2 and ΔPO 2 . We used the Matlab wavelet cross-spectrum toolbox developed by Grinsted et al. [32]. Squared wavelet coherence between the time series of each RGE metric and ΔCBFv was separately plotted with x-axis as time and y-axis as a scale which has been converted to its equivalent Fourier period. An example of squared wavelet coherence between bER and ΔCBFv in right MCA from a representative subject during spontaneous breathing is shown in S3B 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 [32]. An arrow indicates the phase angle between the two time series, with bER leading ΔCBFv, at particular samples of the time-frequency plane: a rightward pointing arrow indicates that the time series are in phase, or positively correlated (ϕ= 0); a leftward pointing arrow indicates anti-correlation (ϕ = π), 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 the Monte Carlo method, and the 5% significance level against red noise is shown as a 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 exceeded 95% significance level, normalized by the maximum possible coherence outside the cone of influence at that particular scale (S3C Fig). It is interpreted in a similar way as the coherence in the transfer function analysis which has been used in cerebral autoregulation study [34].
Interpreting the time-averaged coherence irrespective of phase lag raised the question that the coherence with positive correlation between two time series (at phase lag of 0±π/2) and the coherence with negative correlation (at phase lag of π±π/2) might be mixed together. Therefore, mean time-averaged coherence at the phase lags of 0±π/2 and π±π/2 were separately averaged across all the subjects who participated in the TCD sessions to explore the Fourier periods/frequency bandwidths that oscillations of ΔCBFv were in synchrony with the time series of each RGE metric (bER, ΔPCO 2 and ΔPO 2 ) when they were at rest.

Part 2
Preprocessing of BOLD-fMRI data. All the BOLD-fMRI data were imported into the software Analysis of Functional NeuroImage (AFNI) [35] (National Institute of Mental Health, http://afni.nimh.nih.gov) for time-shift correction, motion correction, normalization and detrending. Details of preprocessing with AFNI are given as follows. 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, motioncorrected and co-registered to the first image of the functional dataset using three-dimensional volume registration. Components of motion were removed from the co-registered functional dataset using orthogonal projection. The clean functional dataset 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 [36,37] (MGH/MIT/HMS Athinoula A. Martinos Center for Biomedical Imaging, Boston, http://surfer.nmr.mgh. harvard.edu) were excluded from the following analyses of functional images. Individual subject brain volumes with time series of percent BOLD signal changes (ΔBOLD) were derived.
Linear regression for the association between ΔBOLD and RGE metrics in individual subject. Linear regression analysis was used to evaluate the association between ΔBOLD, bER, ΔPO 2 and ΔPCO 2 for each subject. For each voxel of the preprocessed brain volume, the time series of ΔBOLD was separately regressed on bER, ΔPO 2 and ΔPCO 2 . Regression coefficient beta (β) value was defined as the percent BOLD signal changes per unit change of the regressor (bER, ΔPO 2 or ΔPCO 2 ). Individual subject brain volumes with β magnitudes were registered onto their own anatomical scans and transformed to the standardized space of Talairach and Tournoux [38]. Monte Carlo simulation was used to correct for multiple comparisons [39]. 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. Based upon a Monte Carlo simulation with 2000 iteration processed with ClustSim program [40], 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 region-of-interest (ROI) analysis for the association between ΔBOLD and RGE metrics. For each subject who participated in MRI scanning, β magnitude values derived by regressing ΔBOLD on bER (β bER ), ΔBOLD on ΔPO 2 (β ΔPO2 ) and ΔBOLD on ΔPCO 2 (β ΔPCO2 ) were separately averaged in each of the 160 brain regions parcellated by the software FreeSurfer. One-sample t-tests were used to analyze regional β bER , β ΔPO2 and β ΔPCO2 in the subject group. False discovery rate was used to correct for multiple comparisons [41,42]. Significant association in group was considered at false discovery rate adjusted p fdr <0.05.
For each brain region, we also calculated the number of voxels with significant ΔBOLD that reached cluster-corrected threshold p<0.05 in the individual subject analysis. The number of voxels with significant ΔBOLD were then normalized to the total number of voxels in each brain region. The percentage of voxels in each brain region with significant ΔBOLD, namely voxelβ bER , voxelβ ΔPO2 and voxelβ ΔPCO2 , for RGE metrics bER, ΔPO 2 and ΔPCO 2 were calculated respectively. Individual subject brain volumes with regional voxelβ due to bER (voxelβ-bER ), ΔPO 2 (voxelβ ΔPO2 ) and ΔPCO 2 (voxelβ ΔPCO2 ) were obtained. One-sample t-tests were used to analyze regional voxelβ bER , voxelβ ΔPO2 and voxelβ ΔPCO2 in the subject group. False discovery rate was used to correct for multiple comparisons [41,42]. Significant association in group was considered at p fdr <0.05.
In the group analysis, we preferred the ROI analysis of 160 brain regions for β bER , β ΔPO2 and β ΔPCO2 instead of the commonly used voxel-by-voxel analysis because the maps of β values were compared side-by-side with the maps of voxelβ bER , voxelβ ΔPO2 and voxelβ ΔPCO2 which require the ROI analysis to evaluate the percentage of voxels in each brain region with significant ΔBOLD.
Group maps of regional ΔBOLD associated with RGE metrics vs. group functional connectivity maps with the seed at left precuneus. Seed-based analysis with the seed at the left precuneus was used to generate the functional connectivity map for an individual subject. For each subject, after preprocessing of BOLD data, a low pass filtering at 0.03Hz was applied onto the brain volume with time series of ΔBOLD. The time series of ΔBOLD in the voxels within left precuneus were averaged. The averaged ΔBOLD time series of left precuneus was correlated with the ΔBOLD time series of each voxel in the brain volume. Pearson's correlation coefficient was used to indicate the strength of correlation and converted to z scores using the Fisher's Z transformation. For the comparison with the regional β and voxelβ maps, the Fisher's z scores in each of the 160 brain regions were averaged for each subject. Individual subject brain volumes with regional Fisher's z scores from all the subjects who participated in MRI sessions were subjected to group analysis using a one-sample t-test. False discovery rate was again used to correct for multiple comparisons [39,41,42]. Significant correlation in group was considered at p fdr <0.05.
Dynamic analysis of coherence between bER and ΔBOLD in brain regions within DMN. To further verify the dynamic coupling between bER and CHF in DMN, we used the WTC method to examine the dynamic changes in coupling between the time series of RGE metrics (bER, ΔPCO 2 , and ΔPO 2 ) and ΔBOLD in three brain regions within DMN. They included inferior parietal lobule (IPL), posterior cingulate (PCC) and precuneus (PCun). The analysis procedures were the same as those described for the dynamic analysis of coherence between ΔCBFv and RGE metrics. The mean time-averaged coherence between each RGE metric and ΔBOLD from each of the three brain regions at the phase lag of 0±π/2 was averaged across all the subjects who participated in MRI sessions. Similarly, the mean time-averaged coherence between each RGE metric and ΔBOLD from each of the three brain regions at the phase lags of π±π/2 were averaged across the same group of the subjects.
Dynamic analysis of coherence between respiratory volume per unit time (RVT) and ΔBOLD in brain regions within DMN. To confirm that the respiratory-variation-related fluctuations reported by Birn et al. [3] had little contribution to CHF within the DMN, we used the same WTC method to examine the dynamic changes in coupling between the time series of RVT (respiration volume per time) and ΔBOLD in three brain regions within DMN (IPL, PCC and PCun). Ten out of 20 subjects had RVT measurements with a respiratory bellow. According to the method described by Birn et al. [3], RVT was computed by the difference between the maximum and minimum bellow positions at the peaks of inspiration and expiration respectively, and this difference was then divided by the period of the respiratory cycle. The mean time-averaged coherence between RVT and ΔBOLD from each of the three brain regions at the phase lag of 0±π/2 were averaged across 10 subjects. Similarly, the mean time-averaged coherence between RVT and ΔBOLD from each of the three brain regions at the phase lag of π±π/2 were averaged across the same group of subjects.
Regional CVR quantification under exogenous CO 2 challenge. Linear regression analysis was used to derive CVR from the time series of ΔBOLD and vasoactive stimulus when the subject was under exogenous CO 2 challenge. We used the time series of P ET CO 2 a regressor in a linear regression analysis. CVR under exogenous CO 2 challenge was defined as the percent BOLD signal changes per mmHg change of P ET CO 2 . Therefore CVR was quantified by the coefficient of regression, i.e. the slope.
For each subject who participated in exogenous CO 2 MRI scanning, CVR values derived from regressing ΔBOLD on P ET CO 2 (CVR CO2-PETCO2 ) were separately averaged in each of the 160 brain regions parcellated by the software FreeSurfer. To study the CVR in group, one-sample t-tests were applied onto the brain volumes with regional CVR CO2-PETCO2 . False discovery rate was used to correct for multiple comparisons.

Results
Subject demographics are shown in Table 1. The mean values and standard deviation of the time series of P ET CO 2 , P ET O 2 , ΔPCO 2 , ΔPO 2 and bER measured in the TCD and MRI sessions are summarized in Table 2.

Correlation among bER, ΔPO 2 and ΔPCO 2
The correlations among the RGE metrics (bER, ΔPO 2 and ΔPCO 2 ) in TCD and MRI sessions are shown in Fig 1 and summarized in S1 Table. Correlation coefficients between ΔPO 2 and ΔPCO 2 ranged from 0.7 to 0.9 in TCD sessions and from 0.2 to 0.9 in MRI sessions, demonstrating that ΔPO 2 and ΔPCO 2 are not necessarily redundant in either TCD or MRI sessions. The range of correlation strength between ΔPO 2 and ΔPCO 2 in the MRI sessions was larger than that in TCD sessions. Correlation between bER and ΔPO 2 was stronger than that between bER and ΔPCO 2 .

Part 1
Correlation between ΔCBFv and RGE metrics. The time series of ΔCBFv in the right MCA of a representative subject during spontaneous breathing with and without applying low pass filtering at 0.03Hz are shown in Fig 2A. The time series of bER, ΔPO 2 , ΔPCO 2 , P ET CO 2 and P ET O 2 acquired simultaneously with ΔCBFv for the same subject are shown in Fig 2B. All the time series including CBFv, bER, ΔPO 2 , ΔPCO 2 , P ET CO 2 and P ET O 2 of the representative subject showed prominent oscillations with periods of 0.5 to 2 minutes which is equivalent to the frequency range of 0.008Hz to 0.03Hz. Oscillations of bER, ΔPO 2 , ΔPCO 2 , P ET CO 2 were in phase with CBFv, while the oscillations of P ET O 2 were out of phase. Oscillation amplitude of ΔPO 2 could double that of ΔPCO 2 (Fig 2B), where the standard deviation values indicating variations over the ΔPO 2 time series were larger than those over the ΔPCO 2 time series as shown in Table 2.
With the simple correlation analyses, Fig 2C shows that the correlation of bER with ΔCBFv was the strongest among all the RGE metrics for the representative subject. ΔPO 2 followed relatively closely behind bER, and ΔPCO 2 was the weakest in correlation with ΔCBFv. ΔPCO 2 and P ET CO 2 shared very similar correlation results with ΔCBFv. The correlation coefficients between ΔCBFv and RGE metrics (bER, ΔPO 2 , ΔPCO 2 and P ET CO 2 ) for all 13 subjects who participated in TCD sessions are shown in S2 Table. The similar correlation of bER and ΔPO 2 with ΔCBFv showed that ΔPO 2 plays an essential role in bER correlation with CHF.
Applying a low pass filter of 0.03Hz to the time series of ΔCBFv (Fig 2A, right panel) increased the strength of correlation of CBFv with bER, ΔPO 2 and ΔPCO 2 (Fig 2C, right  panel) but did not change the order of which RGE metric was more correlated with ΔCBFv.
In the group analysis, paired comparisons of Fisher's z scores transformed from Pearson's correlation coefficients between ΔCBFv and RGE metrics for all 13 subjects are shown in Fig  2D and S2 Table. The Fisher's z scores of the correlation between ΔCBFv and bER were used as a reference to compare with those of the correlation between ΔCBFv and other RGE metrics beside bER. ΔCBFv consistently showed a stronger correlation with bER and ΔPO 2 than with ΔPCO 2 and P ET CO 2 . The correlation of ΔCBFv in LMCA with bER was significantly stronger than that with ΔPO 2 , while the difference between the correlation of ΔCBFv in RMCA with bER and that with ΔPO 2 did not reach statistical significance.
The differences revealed in paired comparisons again demonstrated that ΔPO 2 and ΔPCO 2 were not necessarily redundant.
Dynamic coherence between ΔCBFv and RGE metrics. The time-averaged coherence analyzed by WTC at the phase lag of 0±π/2 which indicates a positive correlation, and at the

PLOS ONE
Dynamic brain-body coupling between resting state hemodynamics and O 2 -CO 2 exchange ratio phase lag of π±π/2 which indicates a negative correlation, were plotted for all 13 subjects participated in TCD sessions (Fig 3). At the phase lag of 0±π/2, the mean time-averaged coherence between bER and ΔCBFv and that between ΔPO 2 and ΔCBFv reached a value of 0.5 or above at the frequency range between 0.008 and 0.03Hz, while the mean time-averaged coherence between ΔPCO 2 and ΔCBFv stayed below 0.1 (Fig 3A, top row). At the phase lag of π±π/2, the mean time-averaged coherence of all three RGE metrics with ΔCBFv stayed below 0.1 at the frequency range between 0.008 and 0.03Hz (Fig 3B, top row).

Part 2
Association between ΔBOLD and RGE metrics in different brain regions. Brain maps on the association between ΔBOLD and RGE metrics (bER, ΔPO 2 and ΔPCO 2 ) during spontaneous breathing in a representative subject are shown in S4A Fig. The regional β magnitudes, namely ΔBOLD per unit change of each RGE metric in different brain regions, were mapped for bER (β bER ), ΔPO 2 (β ΔPO2 ) and ΔPCO 2 (β ΔPCO2 ) (S4A Fig). Comparing with the maps of β ΔPO2 or β bER , the map of β ΔPCO2 showed fewer pixels which reached statistical significance (corrected p<0.05).
In addition to the brain maps from the individual subject, group regional β magnitudes were mapped for bER (β bER ), ΔPO 2 (β ΔPO2 ) and ΔPCO 2 (β ΔPCO2 ) for all 20 subjects who participated in MRI sessions (S4B Fig and Fig 4A). ΔBOLD was found to be significantly associated with bER or ΔPO 2 in most of the brain regions ( Fig 4A) while association found between ΔBOLD and ΔPCO 2 in a few brain regions (insula, medial orbitofrontal and temporal) (S4B Fig) did not reach statistical significance after correcting for multiple comparisons (Fig 4A).
In addition to using the β magnitude, brain volumes with the regional percentage of voxels having significant ΔBOLD for the regressors bER (voxelβ bER ), ΔPO 2 (voxelβ ΔPO2 ) and ΔPCO 2 (voxelβ ΔPCO2 ) were analyzed in the same group of 20 subjects. Significant association of ΔBOLD with bER and ΔPO 2 was shown in more than 50% of the voxels in the gray matter

PLOS ONE
Dynamic brain-body coupling between resting state hemodynamics and O 2 -CO 2 exchange ratio regions (Fig 4B), while the percentage of significant voxels was smaller in white matter regions. A smaller percentage of gray and white matter voxels showed a significant association of ΔBOLD with ΔPCO 2 . The maps of paired comparison between voxelβ bER and voxelβ ΔPO2 showed significant difference in the percentage of significant voxels (Fig 4C) in the rostrum and genu of the left corpus callosum. However, the paired comparison between voxelβ bER and voxelβ ΔPCO2 showed a significant difference in the percentage of significant voxels in many areas, including subcortical brain regions, insula and brainstem.
Regional association between bER and ΔBOLD overlapped with many areas of the default mode network. Group statistical parametric maps of regional association between ΔBOLD and RGE metrics in gray matter and brainstem ( Fig 5A) were used to compare with the group connectivity map with the seed at left precuneus (Fig 5B) for the same group of 20 subjects in the MRI study. Brain regions showing significant association of ΔBOLD with bER and ΔPO 2 included precuneus, posterior cingulate, anterior insula, caudate nucleus, superior temporal and inferior parietal regions. These brain regions overlapped with those in DMN reported by Raichle et al. [23] as well as those reported by Yeo et al. [43] and Choi et al. [44]. Additional regions such as posterior insula, putamen and occipital regions found in our group statistical parametric maps were also shown in the findings by Raichle et al. [23].
Group analysis of individual connectivity maps demonstrated significant increased connectivity in brain regions similar to the regions showing a significant association between ΔBOLD

PLOS ONE
Dynamic brain-body coupling between resting state hemodynamics and O 2 -CO 2 exchange ratio and bER (Fig 5). The close match between Fig 5A and 5B from the same BOLD datasets but acquired with independent methods supports that bER is capable of outlining the connectivity among brain regions in the DMN.

PLOS ONE
Dynamic brain-body coupling between resting state hemodynamics and O 2 -CO 2 exchange ratio Dynamic coherence between bER and ΔBOLD in brain regions within DMN. We used WTC analysis to verify the dynamic coupling between bER and BOLD signal fluctuations in DMN regions such as left and right inferior parietal lobule (IPL), left and right posterior cingulate (PCC), and left and right precuneus (PCun) in all 20 subjects participated in MRI sessions. At the phase lag of 0±π/2, we showed that the mean time-averaged coherence between bER and ΔBOLD and that between ΔPO 2 and ΔBOLD reached a value of 0.2 or above in general (Fig 3A, lower panel) at the frequency range between 0.008 and 0.03Hz, while the mean timeaveraged coherence between ΔPCO 2 and ΔBOLD peaked around 0.2 or below (Fig 3A, lower  panel). At the phase lag of π±π/2, the mean time-averaged coherence of all three RGE metrics with ΔBOLD stayed below 0.1 at the frequency range between 0.008 and 0.03Hz (Fig 3B, lower  panel).
Dynamic coherence between RVT and ΔBOLD in brain regions within DMN. To confirm that the RVT had little contribution to CHF within the DMN, we used the WTC method to examine the dynamic coupling between RVT and ΔBOLD extracted from IPL, PCC and PCun in 10 subjects who had their respiratory tracings measured by a respiratory bellow. We showed that the mean time-averaged coherence between RVT and ΔBOLD at the phase lag of 0±π/2 increased between 0.016 and 0.031Hz, and reduced at the frequency below 0.016Hz where the mean time-averaged coherence started to increase at the phase lag of π±π/2 (S5 Fig). The frequency distribution of coherence of CHF with RVT is shown to be very different from that with bER.
Regional CVR maps under exogenous CO 2 challenge and during spontaneous breathing. Under exogenous CO 2 challenge, most of the brain regions showed increased CVR CO2--PETCO2 in the subject group, especially thalamus, insula and putamen (S6 Fig). During spontaneous breathing, increased β bER was found in most of the brain regions, while no significant association indicated by β ΔPCO2 were shown in most of the brain regions. The CVR map during spontaneous breathing indicated by β bER resembled the CVR map under exogenous CO 2 challenge indicated by CVR CO2-PETCO2 .

Discussion
There has been no previous study on the dynamic interactions between low frequency fluctuations of resting state CHF (ΔCBFv and ΔBOLD) and those of all three RGE metrics of bER, ΔPO 2 and ΔPCO 2 . However, results in this study support such interactions even though the physiological mechanisms are still unclear. The time series of ΔPO 2 and ΔPCO 2 were shown to be interdependent but non-redundant. The prominent oscillations of bER, ΔPO 2 and ΔPCO 2 were characterized by periods of 0.5 to 2 minutes (0.008-0.03Hz) and were coherent with CHF in this frequency range. The coherence of CHF with bER or with ΔPO 2 was much stronger than the coherence of CHF with ΔPCO 2 in the frequency range between 0.008Hz and 0.03Hz. bER was shown in some brain regions to be superior to ΔPO 2 in the association with CHF. Brain regions with the strongest bER-CHF coupling overlapped with many areas of DMN. Taken together, our study provides evidence that resting state CHF at low frequency show the strongest association with fluctuations attributed to bER.

Oscillations of ΔPO 2 and ΔPCO 2 are interdependent but not necessarily redundant
The ΔPO 2 and ΔPCO 2 were often considered to be redundant [4,6,16]. In the present study, the average value of the correlation strength of around 0.8 and the relatively large spread of the correlation between ΔPO 2 and ΔPCO 2 demonstrated that while ΔPO 2 and ΔPCO 2 were interdependent they were not redundant to each other (Fig 1). Respiratory data during spontaneous breathing in healthy infants and adults were consistent with ours, showing that P ET O 2 and P ET CO 2 were correlated but not redundant [13,17,18].
The interdependent but non-redundant oscillations of ΔPO 2 and ΔPCO 2 during spontaneous breathing at the normoxic level are also supported by the findings of feedback loops involving the interaction among chemoreceptors, and blood gases demonstrated in many earlier studies on animal models [10,[45][46][47][48]. Most of these studies on the role of PO 2 to stimulate peripheral chemoreceptors at the carotid body had been focused on hypoxia at PO 2 level less than 60 mmHg [10,49] where chemoreceptor activities rose quickly in a hyperbolic fashion. However, Biscoe et al. [45] showed that peripheral chemoreceptor activities were present from normoxia to hyperoxia up to arterial PO 2 level of 190 mmHg and beyond. Lahiri et al. [48] reported that the stimulus thresholds of arterial PO 2 and PCO 2 for peripheral chemoreceptors were largely interdependent under normoxic conditions. 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 . While research on the mechanisms of interaction between peripheral and central chemoreceptors to optimize systemic blood gases is on-going [9][10][11], the chemoreceptor activities reported in these studies highlight the interdependence of PO 2 and PCO 2 which work synergistically to regulate the blood supply to the brain. bER can be one of the appropriate metrics to be used in the study of such synergism between ΔPO 2 and ΔPCO 2 . The superiority of bER over ΔPCO 2 in correlating with CHF is likely to be attributed partly to the ratio format of bER enabling a reduction of ventilatory fluctuations common to ΔPO 2 and ΔPCO 2 , and partly to the physiological role of bER which takes into account the interaction between ΔPO 2 and ΔPCO 2 . Besides, a regression model of ΔCBFv or ΔBOLD vs RGE metrics, including both ΔPCO 2 and ΔPO 2 as regressors would suffer from collinearity due to the correlation between ΔPCO 2 and ΔPO 2 .
Interestingly, the spread of correlation strengths between ΔPO 2 and ΔPCO 2 as well as those between bER and ΔPO 2 were different between TCD and MRI sessions (Fig 1). Correlation strengths might change in response to different factors which could include the postures of sitting up vs lying down, the noises, or the level of stress from the surroundings. Participants were sitting in an open and quiet environment in TCD sessions while they were lying down in supine position in a noisy environment of the MRI scanner bore. A change from the supine posture to sitting upright was reported to be associated with a redistribution of both blood flow and ventilation in the lungs, which affected the arterial PO 2 [50][51][52][53]. While it is interesting to observe the differences in the correlation strength among RGE metrics between TCD and MRI sessions, the details of such mechanisms are outside the scope of our current study.

Dynamic coupling between CHF and RGE metrics of bER, ΔPO 2 and ΔPCO 2
Correlation and coherence analyses showed the dynamic coupling between CHF (ΔCBFv and ΔBOLD) and RGE metrics of bER, ΔPO 2 and ΔPCO 2 in both TCD (Figs 2 and 3) and MRI sessions (Figs 3 and 4). Had ΔPO 2 and ΔPCO 2 been redundant, their correlation and coherence with CHF would be expected to be the same, and there should be weak correlation and coherence of CHF with bER. In contrast, comparing with ΔPO 2 and bER, we found that ΔPCO 2 had moderate to weak correlation with ΔCBFv (Fig 2). The map of β ΔPCO2 for a single individual showed fewer pixels with a significant association between ΔBOLD and ΔPCO 2 (S4A Fig). Group analysis of regional β ΔPCO2 showed however almost a "null" map as only a few brain regions (insula, medial orbitofrontal and temporal areas) showed association between ΔPCO 2 and ΔBOLD (S4B Fig) but such association did not reach statistical significance after correcting for multiple comparisons (Fig 4). Our group averaged result was consistent with the findings in the study by Golestani et al. [16], where large inter-individual variability was demonstrated in the brain maps of regressing BOLD signal changes on P ET CO 2 during spontaneous breathing.
The interaction between ΔPCO 2 and CHF is accompanied by an interaction between ΔPO 2 and CHF. However, the mechanisms behind the superiority of ΔPO 2 over ΔPCO 2 in association with CHF, especially at the frequency range of 0.008-0.03Hz, remain unclear. A possible factor may be related to the difference in the O 2 and CO 2 dissociation curves, which are associated with the magnitude of the PO 2 and PCO 2 variations in time [13]. During spontaneous breathing, the PO 2 varies from approximately 150mmHg in inspired air to 100mmHg in expired air at the top horizontal asymptote of the sigmoidal O 2 dissociation curve, where it is associated with only a small change in oxygen saturation. On the contrary, the CO 2 dissociation curve is more linear. With the Haldane effect, O 2 affects the affinity of hemoglobin for CO 2 . Given the small fluctuations of oxygen saturation, a rather significant change in PO 2 can, therefore, be associated with a minor change in PCO 2 . The larger oscillation amplitude of ΔPO 2 in comparison with that of ΔPCO 2 was observed by our team (Fig 2 and Table 2) and by Lenfant et al. [13]. An additional factor that was suggested by Lenfant et al. is the variation in the distribution of the ventilation-to-perfusion ratio (V/Q) in the lungs. Lenfant et al. reported that the alveolar-arterial difference in PO 2 decreased as V/Q increased, while that in PCO 2 remained nearly constant and independent of the changes in V/Q. The temporal oscillation of the V/Q distribution could then cause a more considerable change in PO 2 than in PCO 2 . That also raises an intriguing idea for future research to examine any possible relationship between CHF and fluctuations in V/Q. The superiority of ΔPO 2 over ΔPCO 2 in association with CHF suggests a potentially important implication that CHF is more sensitive to the demand of the supply of O 2 than of the removal of CO 2 . In addition to our studying of healthy subjects, exploring the interaction between CHF and RGE metrics in patients with known disorders may increase the understanding of the stronger coupling of CHF with bER and ΔPO 2 than with ΔPCO 2 .

Brain regions with the strongest bER-CHF coupling overlapped with many regions of default mode network
We used the resting state connectivity method with the seed at the left precuneus to outline the areas within DMN (Fig 5B). We found that brain regions with the strongest association between bER and ΔBOLD ( Fig 5A) overlapped with many areas of DMN. Such strong association may be attributed to DMN showing higher metabolic and hemodynamic activities at rest than other parts of the cortex [23]. Fluctuations of systemic gases can be under the influences of ambient (exogenous) gases and/or systemic (endogenous) gases due to metabolism as well as feedbacks between chemoreceptors. Since the brain is one of the major shareholders of systemic blood flow [54] and metabolism [55,56] at rest, any mechanisms that help to preserve cerebral metabolism and CBF within the normal range in the homeostatic regulatory process [12] may leave a signature in the CHF.
In our resting state experiments, CHF were measured on participating individuals during spontaneous breathing. No cognitive, motor, sensory, visual and gas challenges were applied to the participants. The internal environment of the body was steady, and major body organs had constant autonomic communication with the brain back and forth. At rest, background cerebral neurovascular coupling activities are characterized by DMN [23], and RGE metrics can be related to basal activities of the whole body. Hence it would be reasonable for RGE metrics to couple better with DMN than with other brain networks. bER which showed the strongest correlation with CHF is therefore expected to interact more with DMN. The bER-CHF coupling in DMN found in our study is consistent with the findings in previous studies that CBF and DMN connectivity were altered in individuals with disrupted systemic metabolism [57][58][59][60][61][62][63][64].
Physiological processes that may be associated with the coherence between CHF and bER at the low frequency range of 0.008-0.03Hz The CHF coherence with RGE metrics at the 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 (0.008-0.03Hz) have been reported to be related to autoregulation of microvasculature, spontaneous rhythmic oscillations in intracranial pressure (ICP), and intrinsic brainstem rhythm that leads to cerebral blood volume modulation [65,66]. B-waves may also independently reflect the neurovascular coupling process by altering vascular diameters to ensure the delivery of O 2 and other circulating metabolites through the contractile properties of pericytes or vascular smooth muscle cells [67][68][69][70]. Since oscillations of B-waves are not always associated with the change of arterial CO 2 [66], it would be interesting to pursue further the relationship between B-waves and all three of our RGE metrics as well as resting state CHF.
The coherence of CHF with RGE metrics at the frequency of 0.008-0.03Hz may also be related to a proposed mechanism to remove cerebral metabolic waste. In the glymphatic model [71], cerebral metabolic waste was hypothesized to be cleared from the brain via the perivascular space by vascular pulsation. The actual type of vascular pulsation remains to be determined. To describe possible mechanisms for glymphatic convection, Kiviniemi et al. [72] identified a very slow pulsation found in cerebrospinal fluid (0.001-0.023Hz), which is again in the frequency range similar to the oscillations of the RGE metrics.
In the peripheral circulation, laser Doppler measurements showed that the endothelial activity oscillated at the frequency of 0.0096-0.021Hz [73,74], which is within the same low frequency range of CHF and fluctuations of RGE metrics with a particular interest on bER based on our results. Is there any direct relationship between oscillatory cycles of peripheral endothelial activity and CHF? Considering the model of diving reflex [75] where an increase in CBF is associated with peripheral vasoconstriction, is it possible that the decrease in peripheral blood flow is mirrored by an increase in CBF in a homeostatic process during spontaneous breathing? Future research can be directed to explore the association between bER, CHF and peripheral (skin) blood flow oscillations at the low frequency range of 0.008-0.03 Hz.
Separate from the B-waves and glymphatic model, we do not consider heart rate (HR) or heart rate variability (HRV) to be a major contributing factor for the role of RGE metrics in CHF. Even though HR and HRV had been reported to contribute to low frequency fluctuations, the peaks of coherence between ΔBOLD and HR/HRV were at the frequency of 0.05Hz or above (period~30-42 seconds) [5,76] which are different from our findings that RGE metrics were coherent with CHF between 0.008 and 0.03Hz.

Ventilatory volume fluctuations are not the primary origin of our finding on the association between CHF and the fluctuations of RGE metrics
We should point out that ventilatory volume fluctuations are unlikely to be the cause of our observed interaction between CHF and RGE metrics. Respiratory variability, namely changes in respiration volume per time or RVT that had been shown to have a strong correlation with the change in respiratory volume [14] had also been discussed as a possible source for DMN activities [3]. However, a number of considerations showed that RVT did not provide the answers to our findings on the differences observed between ΔPO 2 and ΔPCO 2 in their interaction with CHF. First, RVT would be expected to have the same effect on the time courses of ΔPO 2 and ΔPCO 2 . Second, bER, being a ratio of ΔPO 2 /ΔPCO 2 , reduces the contribution of ventilatory volume fluctuations to the interaction between bER and CHF. Previous studies showed that the relationship between the time courses of P ET CO 2 and RVT was unclear [14,16]. There was only a mild coherence between P ET CO 2 and RVT [14]. The authors in these studies showed that changes in RVT had a weaker correlation with BOLD signal changes when compared with P ET CO 2 [16,77]. Such results are not surprising as P ET CO 2 is supposed to play a more direct role than ventilation to modulate CHF [3,4]. Besides, the time course of RVT is less than ideal as the numerator of RVT acquired with a respiratory bellow indicates the chest excursion, while respiratory volume is the volume of air breathing in or out. Our coherence analysis showed that ΔBOLD at DMN regions were less coherent with RVT than with bER or ΔPO 2 (S5 Fig). RVT showed coherence with ΔBOLD at the phase lag of 0±π/2 between 0.016 and 0.031Hz but the coherence increased slightly between 0.008 and 0.016Hz at the phase lag of π±π/2. The coherence findings of RVT were also different from those of ΔPCO 2 . The frequency bandwidths and the phase differences seen in our coherence of RVT with ΔBOLD are consistent with those reported by Van den Aardweg [14] even though his team was measuring respiratory volume from a facemask. While coherence frequency bandwidths at different phase lags may be referring to the responses from different chemoreceptors (e.g. central vs peripheral), the relationship between the time courses of P ET CO 2 and RVT is still unclear.

Technical issues that have been addressed
In this study, we used TCD to acquire CHF in addition to using fMRI. One advantage of TCD is its high temporal resolution with less concern on the aliasing effect of high frequency hemodynamic signal. Potential technical confounds like gas sampling at two gas analyzers and technical delay on the time series were also addressed in the data acquisition protocol. In the gas sampling circuit, the same gas sample volume was used in both CO 2 and O 2 analyzers at the same gas sampling flow rate. The technical delay on the time series of gas measurements due to transit time of respiratory gases and the response of the equipment were corrected in the preprocessing of the physiological signals.
In the preprocessing of fMRI data, we did not apply RETROICOR [78] to reduce ventilatory and cardiac 'interference' in the BOLD signals because the heart rate was reported to contribute to BOLD signals mainly at 0.3Hz or above [76] which was above the frequency bandwidths of interest (below 0.03Hz) in this study. Our findings of the correlation between ΔBOLD and RGE metrics in healthy subjects were supported by the CBFv data which were acquired with TCD at a sampling rate of 100Hz on the same group of subjects. There is no concern of aliasing of high frequency cardiac signal into our CBFv signal at low frequency range below 0.03Hz. In Fig 2, we also demonstrated that the preprocessing step for CBFv signals like data smoothing did not change the correlation ranking of ΔCBFv with bER, ΔPO 2 and ΔPCO 2 .
In our fMRI regression analysis, we did not convolve the RGE metrics with a specific respiratory response function (RRF) which was commonly used to convolve with RVT or P ET CO 2 to improve its correlation with ΔBOLD [3,4,6,16]. The reason is that RRF is unnecessary when our primary objective is comparing the association of different RGE metrics (bER, ΔPO 2 and ΔPCO 2 ) with ΔBOLD. If the application of RRF improves the correlation of ΔPCO 2 with ΔBOLD, it will also improve the correlation of ΔPO 2 or bER with ΔBOLD. There is no justification for applying RRF to only ΔPCO 2 and not to ΔPO 2 . Applying the same RRF to all of the RGE metrics should have little effect on how one RGE metric ranks against other RGE metrics in its association with ΔBOLD. Another approach to state that convolution of RGE metrics with RRF is not applicable in our case because our research question is 'Which RGE metric contributes more to ΔBOLD?' whereas studies that took up RRF [3,4,6,16] focused on a separate question of 'Which brain regions have stronger ΔBOLD responses to RVT or P ET CO 2 ?' With a different research question, it was reasonable for those previous studies to focus on convolving the RVT or P ET CO 2 time series with RRF which identified the portion of the ΔBOLD that was not directly accounted for by the time series of RVT or P ET CO 2 .

Technical issues that need to be addressed in the future
Even though ΔPO 2 is related to the partial pressure of O 2 utilized in the whole body, we cannot simply substitute ΔPO 2 and ΔPCO 2 for O 2 uptake (VO 2 ) and CO 2 release (VCO 2 ). O 2 uptake and CO 2 release require extra elements like respiratory minute volume (sometimes normalized for body weight) and are typically evaluated as a steady state of RGE data obtained by averaging the breath-by-breath signal over minutes. A measurement of breath-by-breath VO 2 and VCO 2 inside the MRI scanner in addition to the collection of ΔPO 2 and ΔPCO 2 was not included in our study. Given the physical constraint of the MRI settings, we collected ΔPO 2 and ΔPCO 2 using a nasal tubing. We will leave the possibility of measuring VO 2 and VCO 2 during MRI to future neuroimaging research. However, under specific conditions, some information extracted from ΔPO 2 and ΔPCO 2 can be related to that from VO 2 and VCO 2 . As mentioned above, bER has the explicit influence on suppressing the ventilatory volume when we are taking the ratio of ΔPO 2 and ΔPCO 2 . We, therefore, consider the possibility of bER or equivalently ΔPO 2 /ΔPCO 2 being a surrogate for breath-by-breath VO 2 /VCO 2 at rest.
As the vasodilatory role of CO 2 to modulate CBF is an entrenched belief, it is natural to gravitate towards the explanation that association between RGE metrics and CHF must be resulting from "the effect of CO 2 which reflects in O 2 ". However, the model of CO 2 being the sole agent to interact with CHF remains at variance with our findings of non-redundancy between ΔPO 2 and ΔPCO 2 in this study. Therefore, further investigations on the roles of bER, ΔPO 2 and ΔPCO 2 in their association with CHF would be warranted.

Potential application of bER in the evaluation of CVR
While bER was prominently coupled with CHF within DMN, bER was also coupled, albeit to a less extent, with CHF in the rest of voxels throughout the whole brain. Therefore bER can be used instead of ΔPCO 2 as a regressor to evaluate CVR to spontaneous breathing (S6 Fig). The question is whether the expected small perturbations provided by RGE metrics during spontaneous breathing could have any clinical utility. The previous study on the successful identification of local vascular deficits of moyamoya patients [79] had already been reported for CVR to P ET CO 2 obtained in spontaneous breathing. In our recent study, we found that bER is more robust than P ET CO 2 in characterizing the cerebral hemodynamic responses to brief breath hold challenge in healthy adults [15]. Future studies will clarify the sensitivity and benefits of imaging local vascular deficits using bER instead of P ET CO 2 from spontaneous breathing.
Imaging of bER-CHF coupling may open up opportunities for potential clinical diagnosis and intervention. The clinical potential of RER at rest has been reported for decades. RER, the reciprocal of the steady state mean value of bER, has been reported to change significantly with age [80] and with diseases attributed to respiratory [81], liver [63], cardiac [82] and neuronal [83] dysfunctions. Instead of focusing on the steady state results, our manuscript centers more on the dynamic fluctuations of bER, as a separate metric of interest. The bER-CHF coupling may offer an alternative approach to map brain regions interacting with systemic homeostatic processes besides resting state connectivity analysis for the connections among brain regions [2] and task-induced negative BOLD responses [23]. Since abnormal DMN can be an indicator of neuronal disorders [84,85], future studies on different patient populations may help clarify the meaning of abnormal bER-CHF coupling. Our brain-body coupling framework may be applied in the investigation of brain responses to voluntary respiratory techniques like meditation, yoga breathing or tai-chi [86][87][88][89][90][91][92][93] to manipulate RGE for intervention purpose. Instrument assisted respiratory techniques like the Bi-level Positive Airway Pressure (BiPAP) [94][95][96] that applies to sleep apnea or Amyotrophic Lateral Sclerosis (ALS) may also be considered.

Conclusion
While the physiological mechanisms are not completely known for the interaction between CHF and RGE metrics of bER, ΔPO 2 and ΔPCO 2 , our findings provide evidence that fluctuations of RGE metrics are associated with resting state CHF at a low frequency between 0.008Hz and 0.03 Hz. Both bER and ΔPO 2 are more superior to ΔPCO 2 in association with CHF while CHF could correlate more strongly with bER than with ΔPO 2 in some brain regions. Brain regions with the strongest bER-CHF coupling overlap with many areas of DMN. In addition to offering a physiological model to characterize the contribution of gas exchange elements in the low frequency resting state fluctuations, our findings provide new directions to the study of brain-body interaction. Assessing the validity of the popular CO 2only hypothesis as well as other alternative hypotheses to explain our results of bER-CHF coupling remains an on-going project.

S1 Fig. Definition of end inspiration and end expiration on the time series of RGE metrics.
A segment of 90-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 during spontaneous breathing in TCD session. Open circles represent end expiration while closed circles represent end inspiration. 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. (TIF) 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, with bER leading ΔCBFv, 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) Time-averaged coherence 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. The mean time-averaged coherence between time series of RVT and ΔBOLD at the phase lags of 0±π/2 and π±π/2 (thick color lines) in the brain regions of the inferior parietal lobule (IPL), posterior cingulate (PCC) and precuneus (PCun) of the left brain (left panel) and of the right brain (right panel) (n = 10). Color shaded areas represent standard error of the mean. Coherence between two time series at the phase lag of 0±π/2 indicates positive correlation, while negative correlation is represented by the coherence at the phase lag of π±π/2. (TIF) S6 Fig. Regional CVR maps under exogenous CO2 challenge and regional association between ΔBOLD and RGE metrics during spontaneous breathing. The CVR map during spontaneous breathing indicated by β bER changes resembled the CVR map under exogenous CO 2 challenge indicated by CVR CO2-PETCO2 . (TIF) S1 Table. Correlation among RGE metrics. Strength of correlation indicated by Pearson's correlation coefficients among respiratory metrics including bER, ΔPO 2 and ΔPCO 2 in all subjects who participated in TCD sessions (n = 13), and those who participated in MRI sessions (n = 20). (DOCX) S2 Table. 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, ΔPO 2 , ΔPCO 2 and P ET CO 2 (n = 13). 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's Z scores transformed from Pearson's correlation coefficients in groups. Numbers in brackets next to mean Fisher's Z scores indicate p values in the paired comparisons. (DOCX)