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.

Physiological changes including PCO₂, PO₂, 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₂ and PO₂ 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₂ and PO₂ 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₂ and PO₂ 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₂ and O₂ 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.

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 gradient-echo 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₂ and PO₂ 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₂, PO₂, ECG and respiration were measured simultaneously with MRI acquisition. All the physiological measurements, including those of PCO₂, PO₂ 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₂ 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₂ challenge. Given that there is significant inter-individual variance in resting P_ETCO₂ [25], resting P_ETCO₂ was assessed in subjects via calibrated capnograph before the exogenous CO₂ challenge. The subject wore nose-clip and breathed through a mouth-piece on an MRI-compatible circuit designed to maintain the P_ETCO₂ within ± 1–2 mmHg of target P_ETCO₂ [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_ETCO₂). The CO₂ 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_ETCO₂. 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₂ hypercapnic challenge lasted 10 minutes.

When the subject had exogenous CO₂ challenge in MRI session, BOLD-fMRI images were acquired using the same EPI sequence for resting state. The PCO₂ and PO₂ 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.

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₂ and PO₂ was corrected by cross-correlating the time series of PCO₂ and PO₂ 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₂ and PCO₂ (S1 Fig). The breath-by-breath P_ETCO₂ and P_ETO₂ were extracted at the end expiration of PCO₂ and PO₂ time series respectively. ΔPO₂ is defined as (inspired PO₂ –expired PO₂) and ΔPCO₂ is defined as (expired PCO₂ –inspired PCO₂). Breath-by-breath O₂-CO₂ exchange ratio (bER) is defined as the ratio of ΔPO₂ to ΔPCO₂ 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₂ to ΔPO₂ over minutes. The product of ΔPO₂ and ΔPCO₂ was not used to evaluate the interaction with CHF because the effects of fluctuations due to ventilation would be exacerbated in ΔPO₂×ΔPCO₂ (S2 Fig).

Simple correlation analyses were applied to the time series of RGE metrics (bER, ΔPCO₂, and ΔPO₂) in pairs. The correlation was considered significant at p<0.05.

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 one-third 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₂ and ΔPO₂ 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₂, and ΔPO₂, as well as the phase lag of ΔCBFv to bER, ΔPCO₂ and ΔPO₂. 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₂ and ΔPO₂) when they were at rest.

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, motion-corrected 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₂ and ΔPCO₂ for each subject. For each voxel of the preprocessed brain volume, the time series of ΔBOLD was separately regressed on bER, ΔPO₂ and ΔPCO₂. Regression coefficient beta (β) value was defined as the percent BOLD signal changes per unit change of the regressor (bER, ΔPO₂ or ΔPCO₂). 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³ 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₂ (β_ΔPO2) and ΔBOLD on ΔPCO₂ (β_Δ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₂ and ΔPCO₂ were calculated respectively. Individual subject brain volumes with regional voxelβ due to bER (voxelβ_bER), ΔPO₂ (voxelβ_ΔPO2) and ΔPCO₂ (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₂, and ΔPO₂) 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₂ 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₂ challenge. We used the time series of P_ETCO₂ a regressor in a linear regression analysis. CVR under exogenous CO₂ challenge was defined as the percent BOLD signal changes per mmHg change of P_ETCO₂. Therefore CVR was quantified by the coefficient of regression, i.e. the slope.

For each subject who participated in exogenous CO₂ MRI scanning, CVR values derived from regressing ΔBOLD on P_ETCO₂ (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.

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₂, ΔPCO₂, P_ETCO₂ and P_ETO₂ acquired simultaneously with ΔCBFv for the same subject are shown in Fig 2B. All the time series including CBFv, bER, ΔPO₂, ΔPCO₂, P_ETCO₂ and P_ETO₂ 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₂, ΔPCO₂, P_ETCO₂ were in phase with CBFv, while the oscillations of P_ETO₂ were out of phase. Oscillation amplitude of ΔPO₂ could double that of ΔPCO₂ (Fig 2B), where the standard deviation values indicating variations over the ΔPO₂ time series were larger than those over the ΔPCO₂ time series as shown in Table 2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Correlation between the time series of ΔCBFv and RGE metrics during spontaneous breathing.

(A) Time series of ΔCBFv in right MCA of a representative subject at rest, without filtering and after applying low pass filter at 0.03Hz. (B) Time series of RGE metrics of bER, ΔPO₂, ΔPCO₂, P_ETO₂ and P_ETCO₂ acquired simultaneously with ΔCBFv of the same subject in (A). (C) Correlation between ΔCBFv without filtering and the RGE metrics (left), and that between ΔCBFv with low pass filtering at 0.03Hz and RGE metrics (right) of the same subject in (A). (D) For all 13 subjects who participated in TCD sessions, paired comparisons of Fisher's z scores between ΔCBFv in left MCA and bER, and those between ΔCBFv in left MCA and other RGE metrics besides bER (left), and the paired comparisons of Fisher's z scores between ΔCBFv in right MCA and bER, and those between ΔCBFv in right MCA and other RGE metrics besides bER (right). Each gray circle represents Fisher's z score from the correlation analysis of ΔCBFv and the parameter shown on the 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.

https://doi.org/10.1371/journal.pone.0238946.g002

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₂ followed relatively closely behind bER, and ΔPCO₂ was the weakest in correlation with ΔCBFv. ΔPCO₂ and P_ETCO₂ shared very similar correlation results with ΔCBFv. The correlation coefficients between ΔCBFv and RGE metrics (bER, ΔPO₂, ΔPCO₂ and P_ETCO₂) for all 13 subjects who participated in TCD sessions are shown in S2 Table. The similar correlation of bER and ΔPO₂ with ΔCBFv showed that ΔPO₂ 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₂ and ΔPCO₂ (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₂ than with ΔPCO₂ and P_ETCO₂. The correlation of ΔCBFv in LMCA with bER was significantly stronger than that with ΔPO₂, while the difference between the correlation of ΔCBFv in RMCA with bER and that with ΔPO₂ did not reach statistical significance.

The differences revealed in paired comparisons again demonstrated that ΔPO₂ and ΔPCO₂ 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 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₂ 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₂ 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).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Time-averaged coherence between time series of RGE metrics and CHF.

The mean time-averaged coherence in the frequency bandwidths from 0.008 to 0.25Hz (A) at the phase lag of 0±π/2, and (B) at the phase lag of π±π/2, were plotted (thick color lines). Color shaded areas represent standard error of the mean. Coherence between two time series at the phase lag of 0±π/2 indicates a positive correlation, while a negative correlation is represented by the coherence at the phase lag of π±π/2. Top panel shows the coherence of RGE metrics with ΔCBFv in LMCA and RMCA in TCD sessions (n = 13). The lower panel shows the coherence between RGE metrics and ΔBOLD in the inferior parietal lobule (IPL), posterior cingulate (PCC) and precuneus (PCun) within DMN in MRI sessions (n = 20).

https://doi.org/10.1371/journal.pone.0238946.g003

Association between ΔBOLD and RGE metrics in different brain regions.

Brain maps on the association between ΔBOLD and RGE metrics (bER, ΔPO₂ and ΔPCO₂) 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₂ (β_ΔPO2) and ΔPCO₂ (β_Δ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₂ (β_ΔPO2) and ΔPCO₂ (β_Δ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₂ in most of the brain regions (Fig 4A) while association found between ΔBOLD and ΔPCO₂ in a few brain regions (insula, medial orbitofrontal and temporal) (S4B Fig) did not reach statistical significance after correcting for multiple comparisons (Fig 4A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Regional association between ΔBOLD and RGE metrics in the MRI sessions.

(A) Group maps showing the regional ΔBOLD per unit change of bER, ΔPO₂ and ΔPCO₂ (β_bER, β_ΔPO2 and β_ΔPCO2). (B) Group maps showing the regional percentage of voxels with significant ΔBOLD per unit change of bER, ΔPO₂ and ΔPCO₂ (voxel_βbER, voxelβ_ΔPO2 and voxelβ_ΔPCO2). (C) Paired comparisons of voxelβ maps. All the maps had been corrected for p_fdr<0.05.

https://doi.org/10.1371/journal.pone.0238946.g004

In addition to using the β magnitude, brain volumes with the regional percentage of voxels having significant ΔBOLD for the regressors bER (voxelβ_bER), ΔPO₂ (voxelβ_ΔPO2) and ΔPCO₂ (voxelβ_ΔPCO2) were analyzed in the same group of 20 subjects. Significant association of ΔBOLD with bER and ΔPO₂ was shown in more than 50% of the voxels in the gray matter 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₂. 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₂ 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].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Comparison of statistical parametric maps derived from regression and connectivity analyses.

(A) Group statistical parametric maps of the association between ΔBOLD and RGE metrics (bER, ΔPO₂ and ΔPCO₂). The white matter was excluded for the comparison with the group connectivity map. (B) Group connectivity map with the seed at left precuneus. ΔPCO₂ was demonstrated to be much weaker than bER and ΔPO₂ in its association with ΔBOLD in the brain regions of DMN. Comparing the top and bottom rows, similar regions of DMN were outlined with two independent methods of bER-CHF coupling and seed-based connectivity analysis. All the maps had been corrected for p_fdr<0.05.

https://doi.org/10.1371/journal.pone.0238946.g005

Group analysis of individual connectivity maps demonstrated significant increased connectivity in brain regions similar to the regions showing a significant association between ΔBOLD 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.

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₂ 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 time-averaged coherence between ΔPCO₂ 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₂ challenge and during spontaneous breathing.

Under exogenous CO₂ 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₂ challenge indicated by CVR_CO2-PETCO2.