This study conformed to the standards set by the latest revision of the Declaration of Helsinki and was approved by the Institution’s Research Ethics Board. Written informed consent was obtained from all participants. Ten healthy control participants with no history of neurological disorders [age range: 22–70; eight males; age mean (SD) 42.4 (18.68)] were recruited.

All images were acquired using a 3-Tesla GE MRI scanner (Signa HDx, GE Healthcare, Milwaukee, WI, United States), using an eight-channel phased array head coil. Participants had at least 10 min at rest in the scanner before the fMRI time series acquisition. High-resolution T1-weighted 3D spoiled gradient echo sequences were acquired with the following parameters: TR = 7.88 ms, TE = 3 ms, flip angle = 12^°, 146 slices, voxel size = 0.85 × 0.85 × 1 mm, matrix size = 256 × 256, and field of view = 22 × 22 cm. Two BOLD fMRI echoplanar images were acquired, one for each resting-state condition with the following parameters: TR = 2400 ms, TE = 30 ms, flip angle = 70^°, 250 volumes, 41 slices, isotropic voxel size = 3.5 mm, matrix size = 64 × 64, and field of view = 24 × 24 cm.

The dampening of breath-to-breath P_ETCO₂ and the end-tidal partial pressure of O₂ (P_ETO₂) was achieved using an automated gas blender that applies sequential gas delivery algorithms targeting resting P_ETCO₂ (RespirAct™, Thornhill Research Inc., Canada; Slessarev et al., 2007; Fierstra et al., 2013; Fisher, 2016) thereby dampening physiologically significant breath-to-breath variation in P_ETCO₂. Subjects breathed via a soft plastic mask sealed to their face using transparent dressing film (Tegaderm, 3M, St. Paul MN, United States) during poikilocapnia and isocapnia. For isocapnia, P_ETCO₂ was controlled by targeting each subject’s resting P_ETCO₂; average 38 mmHg (range 32–42 mmHg). P_ETO₂ was targeted at 100 mmHg. During the poikilocapnic resting-state scan, participants breathed room air with no targeting of P_ETCO₂ while P_ETO₂ permitting physiological fluctuations as occurs in fMRI experiments. The sequence of isocapnia and poikilocapnia was randomized. All subjects indicated that they were unable to distinguish between the two resting-state conditions.

Neuroimaging data were preprocessed and analyzed using SPM12 (The Wellcome Department of Cognitive Neurology, London, United Kingdom)¹ running in Matlab v7.14 (The Mathworks Inc., United States). All functional volumes underwent slice-timing correction using sinc interpolation to temporally align the slices within each volume. Each volume was spatially realigned to the first volume acquired using a six-parameter rigid body transformation. Realignment motion parameters were set at 3 mm for translation and 1° for rotation relative to the first volume for exclusion criteria. Head movement was measured in three axes with 1.5 mm limits with all subjects meeting this criterion. High-resolution T₁-weighted anatomical volumes were segmented into GM, WM, and cerebrospinal fluid (CSF), and were normalized to Montreal Neurological Institute (MNI) space using the normalized EPI image in SPM. Functional images were normalized into MNI space and smoothed with a spatial convolution 8 mm full-width half-maximum Gaussian kernel.

Functional connectivity measures were analyzed in Matlab using the CONN-fMRI functional connectivity toolbox v17f (Whitfield-Gabrieli and Nieto-Castanon, 2012).² Data were de-noised following the anatomical CompCor approach (Behzadi et al., 2007). CompCor is advantageous in that it does not require external monitoring of physiological fluctuations (cardiac and respiratory) as compared to other noise correction routines. Anatomical CompCor factors in the signal from WM, ventricles, large vessels, and CSF to accurately model physiological fluctuations in GM areas and uses the five most significant principal components each from the WM and CSF as covariates in a general linear model (GLM) as an estimate for physiological noise (Behzadi et al., 2007). We additionally analyzed the resting-state data using CompCor and with CompCor turned off (i.e., removing WM and CSF covariates). We expected that isocapnia would mostly remove variations due to P_ETCO₂ fluctuations, whereas CompCor would remove other sources of noise such as cardiac pulsations, in order to see how these two methods compared.

Global signal regression was not performed in order to bypass introduction of artifactually negative correlations into resulting connectivity measures (Chai et al., 2012). Instead, a number of confounds were regressed out via CompCor to further remove unwanted BOLD signal artifactual effects that were shown to improve sensitivity, specificity, and validity for subsequent functional connectivity analyses. In addition to the CompCor covariates, realignment motion regression (12 regressors comprised of six motion regressors and six first-order temporal derivatives), and functional outliers that were detected via the ART-based identification of outlier scans for scrubbing, thresholded at 0.9 mm for framewise displacement (Power et al., 2012) were also included. As recommended, band-pass filtering was also performed with a frequency window of 0.008–0.09Hz. This preprocessing step helps increase retest reliability, while reducing the effects of low frequency drift and high frequency noise (Weissenbacher et al., 2009). Linear detrending was performed as another de-noising parameter to remove linear trends within each functional scan.

End-tidal PCO₂ was shifted to the maximum correlation with the average BOLD signal using Matlab. Correlation maps between P_ETCO₂ and the resting BOLD signal were computed for each resting-state condition using AFNI software (National Institutes of Health; Cox, 1996) with quadratic fitting to control for baseline and trend. Resulting images were reregistered to MNI space. After Fisher’s Z transformation of the correlation values, group analysis was computed for each group to compare correlations between the P_ETCO₂ time course and BOLD signal during both resting-state conditions. Threshold-Free Cluster Enhancement (TFCE) was performed on the correlations between the BOLD signal and P_ETCO₂ time course data. TFCE is an optimal method in enhancing cluster like structures that has been shown to provide better sensitivity than other methods throughout an extensive range of test signal shapes, signal-to-noise ratios (SNRs), and has been described in detail previously (Smith and Nichols, 2009). Permutation testing was then applied to the height of the maxima of the resulting statistic image, using the “randomize” permutation-based inference tool (Winkler et al., 2014) in FSL v.5.0.9 (FMRIB Library)³ that allowed for the maintenance of strong control over family-wise error (FWE). The effect of spontaneous breathing (poikilocapnia) and (isocapnia) were tested for multiple comparisons, where any significant regions at p < 0.05 were reported and illustrated (Figure 1).

An ROI-to-ROI and Seed-to-Voxel first-level bivariate correlation, hrf-weighted analysis was conducted using a functional connectivity (weighted GLM) model for measuring resting-state connectivity.

The first exploratory whole-brain ROI-to-ROI functional connectivity analysis was performed to examine an unbiased investigation of two types of ROIs within CONN tool. This consisted of a group of 132 ROIs and the second group of 31 ROIs across eight networks that were commonly assigned across subjects. The 132 ROI group was based on the FSL Harvard-Oxford Atlas including 91 cortical, 15 subcortical ROIs [developed at the Center for Morphometric Analysis (CMA), and distributed with the FMRIB Library in FSL], and 26 cerebellar parcellations from the Automated Anatomical Labeling (AAL) Atlas (Tzourio-Mazoyer et al., 2002). The 31 ROI group was derived across eight networks e.g., default mode network (DMN); medial prefrontal cortex (MPFC); posterior cingulate cortex (PCC); right lateral parietal (RLP); and left lateral parietal (LLP) cortex areas. For each ROI, the mean BOLD signal time series was extracted within the ROI voxels, and bivariate correlation coefficients were computed for each pair of ROIs that were Fisher z-transformed. For second-level analysis, a two-sided within-subject paired t-test for all subjects between-conditions contrast of isocapnia vs. poikilocapnia was applied using the false discovery rate (FDR) correction for multiple comparisons p < 0.05. Only significant ROI-to-ROI results are illustrated and reported.

The second hypothesis driven seed-to-voxel analysis was performed using 10 mm spheres positioned in ROIs that have been associated with respiration and rs-fMRI that included the cerebellum, insular/orbitofrontal/cingulate/precuneus/prefrontal cortex, putamen, caudate, superior temporal/supramarginal gyri (Birn et al., 2006; Chang and Glover, 2009), and visual cortex (Madjar et al., 2012). The mean time series were extracted in each seed region by averaging across all voxels for each participant. Bivariate correlation coefficients were computed between the seed time-course and with every other voxel in the brain. The subsequent whole-brain correlation maps were Fisher z-transformed producing z-value maps of voxel-wise functional connectivity for each seed ROI. For second-level analysis, a two-sided within-subject paired t-test for all subjects between-conditions contrast of isocapnia vs. poikilocapnia was applied for voxel-wise statistics throughout the whole brain at an uncorrected level (p < 0.001) before FWE correction was applied at the cluster level (p < 0.05) for multiple comparisons. Only significant seed-to-voxel results are illustrated and reported.

The SD of the breath-to-breath spontaneous P_ETCO₂ variability was significantly lower in the clamped condition 0.84 (0.46) mmHg compared to the unclamped condition 2.47 (0.96) mmHg [mean of variability (SD)], T(9) = 6.06, p = 0.0002 (Figure 2). The mean (SD) P_ETCO₂ was higher in the isocapnic condition 37.53 (3.04) mmHg compared to the poikilocapnic condition 33.81 (3.3) mmHg, T(9) = 8.38, p < 0.0001. This is very unlikely to be of physiologic significance for the following reason. In the isocapnic condition, the P_ETCO₂ is equal to the PaCO₂ (Ito et al., 2008; Fisher et al., 2016). In healthy people P_ETCO₂ is typically 2–4 mmHg less than PaCO₂ due to a small volume of physiological alveolar deadspace (Ito et al., 2008).

Exploratory whole-brain functional connectivity assessment within a whole-brain parcellated network of 132 ROIs and 42 networks revealed significant differences between the two resting state conditions. The controlled breathing isocapnic condition had enhanced functional connectivity between the right putamen with the left intracalcarine cortex (ICC) and the right putamen with the right supracalcarine cortex (SCC) (Figure 3A; Table 1). The poikilocapnic spontaneous breathing condition had enhanced functional connectivity between the right opercular part of the inferior frontal gyrus (IFG) with the left Cerebellar area 1 and the right opercular part of the IFG with the right Cerebellar area 9 (Figure 3B; Table 1). There were no differences between turning CompCor on or off for ROI-to-ROI results.

In addition, a hypothesis-based seed-to-voxel analysis for contrasts comparing isocapnic to poikilocapnic conditions was calculated using second-level GLM analyses. Reported clusters were thresholded at p < 0.05, FWE-corrected, T maxima = 4.78. Enhanced functional connectivity in the isocapnic condition was found between the left putamen seed and a cluster in the precuneus cortex as well as the right putamen seed with a cluster in the left occipital pole and cluster in the inferior division of the right lateral occipital cortex (Figure 4; Table 2). Enhanced functional connectivity in the poikilocapnic condition was found between the primary visual network (PVN) seed and the left hippocampal formation; the left ICC and right SCC seeds with clusters in the left hippocampus; the left SCC seed with clusters in the left hippocampus/lingual gyrus, and callosal body; the posterior cerebellar network (PCN) seed with a cluster in the right supramarginal gyrus (SG); and the cerebellum region 9 seed with the a cluster in the right superior parietal lobule (SPL) (Figure 4; Table 2).

Correlation maps between the P_ETCO₂ and the resting BOLD revealed stronger correlations in the poikilocapnic resting state condition, where P_ETCO₂ was allowed to vary spontaneously as compared with the isocapnic P_ETCO₂ resting state condition. These regions were found predominantly in cerebellar, frontal, and occipital regions (Figure 5). In addition, statistical testing using TFCE for multiple comparison correction p < 0.05 was applied, where significant regions of correlations between P_ETCO₂ with resting BOLD are shown in red (Figure 1). Spontaneous breathing during the poikilocapnic condition showed significantly more P_ETCO₂ and BOLD correlations in GM cerebellar, frontal, parietal, temporal, and occipital regions (Figure 1A) compared to the isocapnic resting state condition (Figure 1B). There were minor differences when using CompCor compared without CompCor in resting state seed-to-voxel analyses.

To the best of our knowledge, this study is the first to compare functional connectivity between poikilocapnic and isocapnic resting-state conditions using whole-brain ROI-to-ROI and seed-to-voxel analysis. The main findings of this study are that implementing isocapnia reveals neuronal functional connectivity within: (1) the corticostriatal pathway between the right putamen and ICC as well as the SCC; and (2) the left putamen and the precuneus cortex. These results are similar to those of Madjar et al. (2012) who found enhanced connectivity between the PCC and occipital regions when P_ETCO₂ was controlled during the resting-state acquisition. Although Madjar et al. (2012) only focused on the PCC as their seed of interest, our findings reveal other resting-state networks that were induced by respiratory changes in PCO₂.

Frontal, temporal, and occipital GM and cerebellum are characterized by large BOLD signal changes in response to a CO₂ challenge (Rostrup et al., 2000; Wise et al., 2004) and therefore may also show the greatest synchronization of BOLD signals to changes in PaCO₂ and hence to each other during poikilocapnia. Functional connections were found between the cerebellum and the IFG and SPL; the cerebellar network with the SG; the visual cortex with the occipital fusiform gyrus; the PVN with the hippocampal formation; and the ICC/SCC with the hippocampus. That these regions are synchronized by P_ETCO₂ fluctuation is shown by loss of synchronization with isocapnia (Figures 1, 5). These findings are consistent with those of Bright et al. (2020) who reported a significant relationship between the more vascularized visual networks during a working memory task under poikilocapnia that was not observed during isocapnia. This finding also supports the notion of a vascular driver for this functional network system as the negative BOLD signal changes were time-locked to the working memory task, whereas visual activation during the working memory task resulted in positive BOLD signal changes.

In addition, a previous study investigating the effects of speech tasks reported decreased CO₂ pressure in the blood (hypocapnia), mainly accounting for the measured changes in oxygenation and cerebral hemodynamics. Their findings suggest the need to monitor P_ETCO₂ during speech tasks as PaCO₂ is still underestimated and a potential confounder in functional brain studies (Scholkmann et al., 2013).

The percentage of BOLD signal change per increment in P_ETCO₂ is lowest in the striatum and thalamus (Posse et al., 1999; Wise et al., 2004), and the basal ganglia (Kastrup et al., 1999). Increased synchronization of the right putamen with the ICC during isocapnia may therefore be explained by withdrawal of confounding vasoactive stimulation, leaving undistorted neuronal synchronization. Thus, vascular responses to changes in P_ETCO₂ that cease with isocapnia are most likely added noise in resting-state functional connectivity analyses and therefore limit sensitivity.

We compared resting-state analyses with and without CompCor noise correction. We found that the significant clusters from seed-to-voxel analysis were slightly larger in volume using CompCor. These results may be indicative of how isocapnic and CompCor operations remove correlated components. We expected isocapnia to mostly remove variations due to P_ETCO₂ fluctuations, whereas CompCor would remove other sources of noise including, for example, those due to cardiac pulsations. Tellingly, isocapnic data were unaffected by application of CompCor. It may be that isocapnia predominantly removes artifactual correlations between highly vasoactive regions, while CompCor predominantly removes artifacts from different spatial locations, revealing when they are actually operating independently. Such post processing should reduce false positives and true positives, whereas isocapnia should not reduce true positives. With control fluctuations of PaCO₂, CompCor may be used to account for respiratory motion and cardiac pulsations, which persist with isocapnia.

It would seem reasonable to regress out the variation due to breath-to-breath variability in P_ETCO₂ and generate cerebrovascular reactivity data expressed as Δ BOLD/Δ P_ETCO₂. Unfortunately, the general correlation between PaCO₂ and P_ETCO₂ breaks down at the scale of small breath-to-breath P_ETCO₂ fluctuations. A small reduction in tidal volume, for example, may result in an incomplete sampling of alveolar gas at the sensor resulting in an erroneously low P_ETCO₂. However, the same breath also results in reduced alveolar ventilation causing an increase in PaCO₂. This breath-by-breath divergence of PaCO₂ from P_ETCO₂ confounds calculations of CVR from changes in P_ETCO₂, while breathing at rest. Moreover, Golestani and Chen (2020) sought to quantify the difference between rs-fMRI measures during isocapnia, poikilocapnia, and poikilocapnia after P_ETCO₂ correction using a CO₂ response function (HRF_CO2) (Golestani and Chen, 2020). Their findings revealed that clamping P_ETCO₂ during isocapnia did not significantly affect respiratory volume variability and cardiac-rate variation, but significantly reduced functional connectivity as calculated with voxel seeds and reduced inter-subject variability in functional connectivity. In addition, applying the P_ETCO₂ correction during poikilocapnia showed to increase inter-subject variability in functional connectivity (Golestani and Chen, 2020).

The small sample size may have limited our ability to detect additional regions masked by respiratory effects. However, that we have found a consistent difference in synchronization of specific brain territories with and without isocapnia despite the inter-subject variability and our small sample size suggests that the effect size of poikilocapnia is large.

Changes in perfusion pressure can affect cerebral blood flow (CBF) during hypercapnia (Aaslid et al., 1989; Perry et al., 2014). An increase in mean arterial pressure (MAP) can alter CBF in the absence of cerebral micro-vessel vasodilation (Claassen et al., 2007). The importance of accounting for MAP would help solve for the magnitude of cerebral vasodilatory response, which is important for adequately matching cerebral perfusion with neuronal metabolic demand. However, MAP cannot be measured breath-to-breath with a non-invasive blood pressure cuff, which would confound the matching of blood pressure to breath-to-breath changes in P_ETCO₂ and BOLD signal. In our study, we did not measure MAP. A common relationship in previous studies showed elevated MAP responses to hypercapnia in older adults compared to young adults (Claassen et al., 2007; Coverdale et al., 2017; Miller et al., 2018). However, in our studies hypercapnia was not applied and there was no reason to suspect synchrony of blood pressure with respiration. Additionally, there may be other mechanisms hampering rs-fMRI brain activity arising from other non-neuronal sources such as contributions from sympathetic vascular innervation, which may affect the fMRI signal through its connection with CBF regulation (Özbay et al., 2019) that warrant further investigation. A theoretical possibility is that eliminating the breath-to-breath variability of arterial PCO₂ would affect the intrinsic connectivity pattern. However, the actual finding in this paper remains that the pattern seen with poikilocapnia includes false synchrony due to vasoactivity. As such, reports of connectivity data collected under poikilocapnia would need to address this issue.

We found that isocapnia decreases or ceases the synchronization in some brain regions previously considered part of the connectome. This suggests that fluctuating P_ETCO₂, during the acquisition of resting-state functional data generate CO₂-synchronized changes that would be considered false positive results for neuronal based connectivity. This artifact is particularly prominent in regions that are highly vascularized and have high vasoactive responses to PaCO₂ such as the visual, cerebellar, and frontal regions. In those, and possibly other regions, the magnitude of the BOLD signal changes resulting from vasoactive stimulation may overwhelm those due to resting state neurovascular coupling, such as in the putamen and calcarine cortex, resulting in false negative connectivity readings in those regions. The findings indicate that resting state acquisitions would benefit from isocapnia in identifying neuronal resting state functional connectivity.

The important finding in this paper remains that the pattern seen with poikilocapnia includes false synchrony due to vasoactivity. Armed with this, the ability to assess how other conditions or diseases alter DMN activity can be more accurately identified and quantitated. This has the added benefit of reducing the sample size in studies trying to determine how experimental conditions, drugs, or diseases influence activity in the DMN.

JF and DM are among the developers of the RespirAct™ for MRI studies at the University Health Network, part of the University of Toronto. Thornhill Research Inc. (TRI) is a for-profit biomedical manufacturing company that was spun off from UHN. It assembles the RespirAct™, on a non-profit basis to enable MRI research at UHN and around the world. JF receives income for work done for TRI and DM holds a minor equity position in TRI.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.