|Home | About | Journals | Submit | Contact Us | Français|
A significant component of BOLD fMRI physiological noise is caused by variations in the depth and rate of respiration. It has previously been demonstrated that a breath-to-breath metric of respiratory variation (respiratory volume per time; RVT), computed from pneumatic belt measurements of chest expansion, has a strong linear relationship with resting-state BOLD signals across the brain. RVT is believed to capture breathing-induced changes in arterial CO2, which is a cerebral vasodilator; indeed, separate studies have found that spontaneous fluctuations in end-tidal CO2 (PETCO2) are correlated with BOLD signal time series. The present study quantifies the degree to which RVT and PETCO2 measurements relate to one another and explain common aspects of the resting-state BOLD signal. It is found that RVT (particularly when convolved with a particular impulse response, the “respiration response function”) is highly correlated with PETCO2, and that both explain remarkably similar spatial and temporal BOLD signal variance across the brain. In addition, end-tidal O2 is shown to be largely redundant with PETCO2. Finally, the latency at which PETCO2 and respiration belt measures are correlated with the time series of individual voxels is found to vary across the brain and may reveal properties of intrinsic vascular response delays.
Functional magnetic resonance imaging (fMRI) is widely employed to study spatial and temporal patterns of neural activity (Bandettini et al., 1992; Kwong et al., 1992; Ogawa et al., 1992). However, the blood-oxygen level dependent (BOLD) signal measured by fMRI is influenced not only by neuronal activity, but also by the effects of cardiovascular and respiratory processes (Biswal et al., 1993; Dagli et al., 1999; Jezzard et al., 1993; Weisskoff et al., 1993). Such physiological noise can obscure neuronal signals because the noise and signals can overlap in time, be anatomically co-localized, and can have (or alias into) common temporal frequencies. It is therefore important to understand how physiological processes – which may be unrelated to neural activation – influence fMRI time series.
The respiratory process introduces noise into fMRI time series in several ways. First, respiration causes modulation of the static magnetic field due to thoracic and abdominal movement, which, if uncorrected, results in distortion, blurring, or displacement of the image, depending on the k-space trajectory. This type of artifact can be reduced prospectively with a navigator correction (e.g. (Pfeuffer et al., 2002)), or retrospectively by filtering out a subspace of the voxel time series that is time-locked to the respiratory cycle (Glover et al., 2000).
Second, natural fluctuations in the depth and rate of breathing over the course of an fMRI scan may cause substantial BOLD signal fluctuations by directly modulating deoxyhemoglobin concentrations. The respiration volume per time (RVT), which quantifies respiratory variations using a pneumatic belt to measure breathing-related chest expansion, was demonstrated to have significant time-lagged correlations with fMRI time series (Birn et al., 2006). It was subsequently shown that convolving the RVT time series with a particular impulse response (the respiration response function (RRF); Fig. 1) yields an optimal linear mapping of RVT to fMRI time series in both cued-inspiration paradigms and task-free resting state (Birn et al., 2008; Chang et al., 2009). According to this model, the average BOLD signal response to a brief increase in RVT exhibits a small initial increase, followed by a larger undershoot occurring approximately 15–20 s later.
Importantly, the mapping between RVT and BOLD signal changes is significant across widespread regions of gray matter, thus impeding the detection of neural activation across many tasks as well as altering correlations between brain regions (functional connectivity). Correcting for the effects of respiratory variation has been shown to improve the detection of task-activated voxels and reduce false positives in resting-state functional connectivity analyses (Birn et al., 2006; Chang et al., 2009).
It is assumed that respiration modulates the BOLD signal primarily by inducing changes in levels of carbon dioxide (CO2), a potent vasodilator (Birn et al., 2006). This assumption is supported on several grounds: (1) the partial pressure of arterial CO2 (PaCO2) is a function of ventilation, which is defined as the product of the depth and rate of respiration (Berne and Levy, 1993); (2) inhaling elevated concentrations of CO2, as well as holding one’s breath, initiates cerebral vasodilation leading to increased cerebral blood flow (CBF) and corresponding BOLD signal changes (Bandettini and Wong, 1997; Kastrup et al., 1999; Kwong et al., 1995; Li et al., 1999; Nakada et al., 2001; Rostrup et al., 2000; Stillman et al., 1995; Vesely et al., 2001; Wise et al., 2007); (3) spontaneous fluctuations in end-tidal CO2 (PETCO2) have been observed to correlate significantly with resting-state fMRI time series (Wise et al., 2004); (4) brain regions that correlated with PETCO2 in (Wise et al., 2004) appeared to be similar to those correlating with RVT in (Birn et al., 2006).
However, the degree to which RVT and CO2 fluctuations relate to one another, or explain common components of the BOLD signal – and across common brain regions – is unknown. If RVT (perhaps via the RRF transfer function) is primarily a predictor of CO2 concentration, it is possible that a more direct measurement of CO2 in the body, such as PETCO2, predicts physiological BOLD signal changes more accurately than RVT. RVT reflects only chest expansion and is at best an approximation of tidal volume, neglecting factors such as dead space and metabolic CO2 output, and thus appears to be a rather indirect measurement of respiratory processes. On the other hand, RVT may account for variance in the BOLD signal beyond that which can be explained by CO2 levels. For instance, expansion/contraction of the chest during breathing elicits a cascade of processes in which CBF is modulated by intrathoracic pressure levels, heart rate, and corresponding autoregulatory processes (Nakada et al., 2001; Thomason et al., 2005). Therefore, although CO2 is a key intermediate variable, it may explain merely a subset of the involved processes.
In the current study, respiration belt measurements (RVT and RVTRRF, where the latter is defined as RVT convolved with the RRF), PETCO2, and the resting state BOLD signal were directly compared. The major aims of this study were (1) to quantify the temporal relationship between RVT and PETCO2 fluctuations, thereby probing the mechanism by which RVT explains BOLD signal changes, and (2) to quantify the temporal and spatial relationship between RVT(RRF), PETCO2 fluctuations, and BOLD signal changes, thereby determining the relative efficacy of RVT(RRF) and PETCO2 in explaining physiological noise across the brain. From a practical standpoint, if a measurement derived from the respiration belt accounts for a BOLD signal component overlapping with that of PETCO2, there may be little additional benefit to monitoring PETCO2 for purposes of reducing physiological noise. In addition, we aimed to (3) compare PETO2 levels to time series measurements from the respiratory belt, PETCO2, and BOLD signal. Although PETO2 changes are likely to mirror PETCO2 changes at rest (via strong negative correlations) and BOLD reactivity to PETO2 is much smaller than that of PETCO2 (Prisman et al., 2008), the effectiveness of spontaneous, resting PETO2 in explaining BOLD physiological noise is unknown and merits investigation.
Participants included 7 healthy adults (3 females, aged 30.9 ± 15.7 years). All subjects provided written, informed consent, and all protocols were approved by the Stanford Institutional Review Board. Data were collected from one additional subject, but were discarded because the CO2/O2 traces did not yield definitive end-tidal values (maxima for CO2 or minima for O2) across the majority of the scan.
Magnetic resonance imaging was performed on a 3.0T whole-body scanner (Signa 750, GE Healthcare Systems, Milwaukee, WI) using an 8-channel head coil. Head movement was minimized with a bite bar. Thirty oblique axial slices were obtained parallel to the AC-PC with 4-mm slice thickness, 0-mm skip. T2-weighted fast spin echo structural images (TR = 3000 ms, TE = 68 ms, ETL = 12, FOV = 22 cm, matrix 192 × 256) were acquired for anatomical reference. A T2*-sensitive gradient echo spiral in/out pulse sequence (Glover and Lai, 1998; Glover and Law, 2001) was used for functional imaging (TR = 2000 ms, TE = 30 ms, flip angle = 77°, matrix 64 × 64, same slice prescription as the anatomic images). A high-order shimming procedure was used to reduce B0 heterogeneity prior to the functional scans (Kim et al., 2002). Importantly, a frequency navigation correction was employed during reconstruction of each image to eliminate blurring from breathing-induced changes in magnetic field; no bulk mis-registration occurs from off-resonance in spiral imaging (Pfeuffer et al., 2002).
Subjects underwent one 16 min fMRI scan. Ten seconds after the beginning of the scan, subjects were cued to take a deep breath and hold it for 6 s. This was performed in order to define a reference point for aligning the respiration, CO2, and O2 waveforms, as described below. For the remainder of the scan, subjects were instructed to rest quietly with their eyes closed.
During the fMRI scan, subjects wore a nasal cannula through which CO2 and O2 were continuously sampled at 100 ml/min using capnography (Capnomac Ultima, Datex Corp., Helsinki, Finland). The cannula was connected to the capnograph with standard Luer fittings and gas sampling tubing of 7 m length. Respiration was monitored with a pneumatic belt placed on the subject’s abdomen just below the rib cage and connected to an arterial line blood pressure transducer (Abbott Laboratories, Chicago: 4285-05) with 6 m of hose. The transducer signal was amplified by a bio-instrumentation amplifier (CB Sciences Inc., Dover N.H.: ETH-250). Analog signals from the capnograph (CO2 and O2) and the respiration amplifier were recorded using a data logger sampling at 100 Hz (Measurement Computing, Norton, MA). For synchronization with the fMRI data, recordings were initiated by a trigger sent from the scanner. Note that the scanner’s built-in respiration monitoring system incorporates a high pass filter and automatic gain control that can distort the respiratory signal. Comparisons between recordings from the home-made belt and the scanner demonstrated that the readings were largely compatible during the conventional breathing episodes employed here, but the former was nevertheless used for purposes of maximum accuracy. Cardiac events were recorded at 100 Hz using the scanner’s built-in photoplethysmograph placed on a finger of the left hand.
RVT was computed from the respiration waveform as described in (Birn et al., 2008). First, peak detection was used to determine the maximum and minimum for each breath. The time series of maxima and minima (max(t) and min(t)) were then interpolated to the imaging TR. A time series of the breathing period (T(t)) was determined by taking the time difference between successive maxima, assigning the resulting value to the midpoint between the successive maxima, and interpolating the resulting time series to the imaging TR. The final RVT time series was then computed as RVT(t) = (max(t) − min(t))/T(t). In addition, RVT(t) was convolved with the RRF (Birn et al., 2008), forming a separate waveform (RVTRRF(t)). Heart rate was calculated from the photoplethysmograph data as the inverse of the average beat-to-beat interval in a 6 s sliding window.
Because there is a delay between a subject’s breathing cycle and the time at which the expired gases reach the capnograph to be analyzed and recorded, the continuous CO2 and O2 data were first aligned to the raw respiration data by maximizing the cross-correlation (for lags between ±10s) with the raw respiration waveform. Data were visually inspected to ensure that the initial breath hold was properly aligned across the 3 waveforms.
PETCO2 was computed from the raw CO2 data by finding the maximum value at each exhaled breath (i.e. lying between successive maxima of the respiration waveform), and interpolating the values to the imaging TR. Similarly, PETO2 was computed by interpolating between the minimum values of the O2 data for each exhaled breath. End-tidal levels of CO2 and O2 have been shown to closely reflect their respective arterial partial pressures (Robbins et al., 1990).
Unfortunately, the capnograph occasionally performed a self-calibration process lasting 30–40 s, during which no CO2 and O2 data were logged. Therefore, for each subject, the longest contiguous segment of data between periods of self-calibration was selected for analysis. Although one might instead have chosen to interpolate across the missing points, the long duration of missing values would impact the cross-correlation analysis (described below). The first 60s of physiological data were also discarded in order to ignore effects of the initial breath hold. The resulting segment of data selected for each subject was at least 5.9 min long (mean=10.9 min, SD=4.0; Table 1). All subsequent analyses of physiological and fMRI data were restricted to this pre-defined window of time. Finally, a linear trend was removed from each physiological time series.
Functional images were pre-processed using custom C and Matlab routines. Pre-processing included slice-timing correction using sinc interpolation, spatial smoothing with a 3D Gaussian kernel (FWHM=5 mm), and removal of linear and quadratic temporal trends. Spatial normalization to a standard template was not performed as a pre-processing step, to avoid blurring; all computations were done in the original subject-space, and voxels were maintained in their original dimensions (3.4375 mm × 3.4375 mm × 4 mm). Data were subsequently corrected for artifacts time-locked to the cardiac and respiratory cycles using RETROICOR (Glover et al., 2000).
Motion parameters were calculated using methods described in (Friston et al., 1996). Image coregistration was not performed because of possible interpolation errors and unintended smoothing across voxels. Therefore, even though a bite bar was used, it was important to verify that motion was minimal. The maximum peak-to-peak excursion and root mean square (RMS) fluctuation for the 3 translational and 3 rotational motion parameter time series were calculated for the segment of resting state data used in the analysis. Rotations were converted to worst-case translations by multiplying by 65 mm, an average head radius (Thomason and Glover, 2008). Summary statistics are reported as the maximum of these values over the 6 axes of motion.
The relationship between physiological signals was examined using pair-wise cross-correlation, allowing time shifts between −120 s to 120 s in increments of 2 s (the imaging TR). This relatively broad range of time shifts was chosen to encompass both immediate influences between physiological processes (such as the causal impact of respiration on heart rate and PETCO2) as well as longer-range influences (such as chemically-mediated feedback between PETCO2 and respiration that may occur on the order of 60 s (Van den Aardweg and Karemaker, 2002)). Maximum cross-correlation magnitudes (rmax) were transformed to Z-scores using the formula , where N is the number of time points. Z-scores are provided in addition to rmax in order to facilitate a comparison across subjects, since a different number of time points (degrees of freedom) were used for each subject.
Physiological signals were individually cross-correlated with every voxel in the brain. Time shifts ranged from −10 s to 30 s in increments of 0.5 s, which represents the range of latency values previously observed between these signals and the fMRI time series (Birn et al., 2006; Birn et al., 2008; Wise et al., 2004). The maximum cross-correlation magnitude (r max), as well as the time lag maximizing the cross-correlation (latency), were determined for each voxel. In accordance with previously-established models (e.g. (Birn et al., 2008; Wise et al., 2004; Wise et al., 2007)), rmax was defined to be the maximum positive cross-correlation magnitude for both PETCO2 and RVTRRF, whereas for PETO2 and RVT, rmax was defined to be the maximum negative cross-correlation magnitude.
For comparison across subjects, rmax values were converted to Z-scores using the formula above. A threshold of Z>5.3 was chosen for purposes of displaying individual subjects’ spatial maps and for comparing the variance explained by (and latency of) physiological signals in the brain. Since cross -correlation involves computing correlation coefficients over a range of time lags, significance testing of the maximum cross-correlation coefficient cannot be carried out the same way as for a correlation coefficient computed at a single time lag. Thus, to guide the choice of this threshold, a Monte Carlo simulation was performed to estimate a null distribution for the maximum cross-correlation amplitude between a respiratory signal and a BOLD signal time series. To generate signals having frequency/autocorrelation characteristics of a typical respiratory signal, independent instances of uniform white noise were filtered with the amplitude of the Fourier transform derived from Subject 1’s PETCO2 signal. Six thousand such instances of random “respiratory” signals were cross-correlated, over a 40-s range of time lags, with the time series of 6000 randomly-selected voxels (1000 voxels from each of the remaining 6 subjects’ brains; voxels were randomly selected in order to minimize spatial proximity, i.e. dependence). A histogram of the maximum cross-correlation magnitudes (converted to Z-statistics) was formed to estimate a null distribution, under which Z=5.3 was found to correspond to p=0.05. Note that for a (single-lag) correlation, Z=5.3 would correspond pond to p<1e−7, demonstrating that cross-correlation can indeed inflate correlation values. We note, however, that the emphasis of the present study is on comparing the relationships between different respiratory measures with the brain, rather than establishing the absolute significance of any one measure.
In previous work, PETCO2 was convolved with a conventional gamma variate impulse response (with a width of 6 s and mean lag of 6.3 s) before assessing its correlation with the brain (Wise et al., 2004), as this impulse response was shown to well-characterize the relationship between PETCO2 and CBF velocity as measured with transcranial Doppler ultrasound (Panerai et al., 2000). While the cross-correlation analysis employed in the present study inherently accounts for variable time delays, convolution with a gamma variate also introduces smoothing, which may alter correlations with the brain. Thus, the cross-correlation was additionally performed using PETCO2 (and PETO2) convolved with a gamma variate function (Glover, 1999), which will be referred to as PETCO2-GAM (PETO2-GAM).
It is of interest to know whether significant additional variance in the BOLD signal can be explained by modeling PETCO2 and RVTRRF simultaneously, beyond that explained by either signal alone. Thus, for each voxel a design matrix containing 2 regressors (PETCO2 and RVTRRF, each time-shifted by the latency maximizing its respective cross-correlation with the voxel’s time series), was formed. Since unconstrained linear regression can assign negative coefficients (betas) to PETCO2 or RVTRRF in order to minimize the mean squared error of the fit, a nonnegativity constraint on the betas was incorporated into the linear regression in order to preserve the interpretability of the model. The amount of variance explained by this simultaneous model was compared to the variance explained by the maximum positive cross-correlations between PETCO2 and the BOLD signal, and between RVTRRF and the BOLD signal. Significance was assessed for each voxel using an F-test (p<0.001). In addition, to depict brain regions most significantly responsive to respiration at the group level, individual subject T-maps from the simultaneous model (both PETCO2 and RVTRRF) were normalized to a common EPI template and entered into a group-level random-effects analysis using SPM5 (www.fil.ion.ucl.ac.uk/spm/software/spm5/).
Lastly, PETCO2 and RVTRRF latency maps (which depict for each voxel the optimum time lag at which PETCO2 and RVTRRF, respectively, were cross-correlated with that voxel’s time series) were compared. The set of voxels included in the comparison was restricted to those which (1) exceeded a threshold of Z>5.3 and (2) had a latency value that did not lie on the upper or lower bounds of the range of cross-correlation lags. These criteria were enforced because the reliability of the latency estimate suffers when the magnitude of the cross-correlation is low, and because it cannot be determined whether latencies lying at the boundaries of the queried interval are local maxima. Latency maps were compared using Pearson correlation. To gauge whether regional patterns in relative latency across the brain emerge at the group level, an average latency map was constructed from spatially normalized single-subject latency maps. To account for inter-subject differences in absolute timing, the latency map of each subject was mean-centered prior to averaging. One subject (Subject 2) was excluded from the averaging because this subject had a strongly bimodal distribution of latencies (a mixture-of-Gaussians fit yielded 2 peaks at −1.2 ± 4.6s and 20.1 ± 4.0s, respectively) whereas the other subjects had nearly unimodal latency distributions spanning a narrower range.
The mean ± standard deviation of PETCO2 and PETO2 signals for each subject are provided in Table 2. Across subjects, the average PETCO2 value was 39.9±3.7 mmHg, and the average PETO2 value was 14.7±0.9 % of room air.
Motion was verified to be minimal. Across subjects, the average drift across the scan session was <0.7 mm (max = 1.06 mm), and the average RMS jitter was only 0.1 mm (max = 0.16 mm).
The RVT, RVTRRF, PETCO2, and PETO2 signals for one representative subject are illustrated in Fig. 2. The maximum cross-correlation coefficient and corresponding time delay between gas (PETCO2, PETO2) and respiration belt (RVT, RVTRRF) waveforms are shown in Table 3, and the cross-correlation magnitudes are represented graphically in Fig. 3.
PETCO2 had strong negative correlations with RVT ( = −14.3, SD=2.66), for which fluctuations in RVT preceded corresponding fluctuations in PETCO2 by an average of 5.14 s (SD=1.57 s). For all subjects, the cross-correlation function had a single clear extremum (a minimum) within queried range (±120s) of time lags. The direction of causality from RVT to PETCO2, along with the negative relationship, is consistent with the ventilatory influence on CO2 levels (Berne and Levy, 1993; Van den Aardweg and Karemaker, 2002).
The cross-correlation between PETCO2 and RVTRRF was positive ( = 18.3, SD=2.55), having a single clear maximum, and was of significantly greater magnitude than cross-correlations between PETCO2 and RVT (p<0.05, paired t-test). PETCO2 preceded RVTRRF by an average of 11.14 s (SD=1.07 s). RVTRRF represents the modeled BOLD signal output to variations in breathing (as measured with a respiration belt), and its high correlation with the (time-delayed) PETCO2 indicates that the effects accounted for by the 2 signals overlap to a large degree (Fig. 2B). We note that since the main feature of the RRF is a wide, negative undershoot, performing a convolution between RVT and the RRF acts primarily to negate and smooth the RVT signal; thus, the increase in correlation magnitude obtained with RVTRRF (compared to RVT) springs primarily from the degree of smoothing introduced by the RRF. Cross-correlation magnitudes between PETCO2-GAM and RVTRRF were not significantly different than those of PETCO2 and RVTRRF (paired t -test, p>0.05).
Correlations between PETCO2 and PETO2 were strongly negative ( = −22.3, SD=3.94), and shifted by either 0 or 1 lag (1 lag = 2 s); see Table 3. As evident from Fig. 3, the correlations between PETO2 and respiratory belt measurements were generally not as strong as those with PETCO2.
Across subjects, the cross-correlation functions between RVT and heart rate, as well as between PETCO2 and heart rate, were variable and showed no obvious trends (data not shown).
Spatial maps of the voxel-wise percentage of BOLD signal variance accounted for by PETCO2 and RVTRRF are presented in Fig. 4(A and B). Variance was explained in gray matter voxels across all subjects, though the extent and magnitude varied considerably across subjects (Table 4, Fig. 5).
Figure 4 indicates that, for a given subject, the spatial distribution of variance explained by PETCO2 is similar to that of RVTRRF. To confirm this impression, a voxel-by-voxel (spatial) correlation coefficient between the rmax maps of PETCO2 and RVTRRF was calculated (Table 5). A scatterplot of the voxel-wise rmax values for PETCO 2 and RVTRRF is shown for one subject (Subject 4) in Fig. 6. Spatial correlations were high for all subjects, indicating that the relative pattern of variance explained across the brain by the two signals is similar. However, differences in the absolute strength with which PETCO2 and RVTRRF related to the BOLD signal were present, as apparent in Fig. 4 and further indicated by Table 4. In some subjects and voxels, RVTRRF showed a stronger relationship with the BOLD signal than PETCO2, while in others the opposite was true.
Spatial maps of the percentage of variance explained when modeling both PETCO2 and RVTRRF together are displayed in Fig. 4C. For comparison with Table 4, Table 6 quantifies the percentage of voxels in the brain (relative to the subject’s whole-brain mask) for which the fitted linear combination of PETCO2 and RVTRRF was correlated above Z>5.3, as well as the average R2 (mea n ± SD) of the fit. Table 6 also provides the percentage of voxels in the brain for which modeling both PETCO2 and RVTRRF together explained significantly more variance (F-test, p<0.001) than either PETCO2 and RVTRRF alone.
As expected, including both PETCO2 and RVTRRF in the model explained more variance than either signal individually. The additional benefit of including both regressors varied across subjects; for example, including both PETCO2 and RVTRRF was significantly better (at the given threshold) than RVTRRF alone in only 5.2% of voxels for Subject 4, but in 38% of voxels for Subject 1. As reflected by differences between Figs. 4A and 4B, for some subjects (e.g. Subject 1), a greater advantage was obtained over PETCO2 alone than over RVTRRF alone, whereas for other subjects (e.g. Subject 3) the opposite was true.
Consistent with the strong negative correlations between PETCO2 and PETO2 (Table 3), both accounted for similar BOLD signal variance and extent (Table 4). It was also observed that convolving PETCO2 (or PETO2) with the gamma variate impulse response function did not effectively alter the cross-correlation magnitude between PETCO2 (or PETO2) and the brain (averaged across subjects, there was no significant difference in variance explained for 99.7% of voxels in the brain; p<0.05, Hotelling’s t-test for dependent correlations (Hotelling, 1940)). On the other hand, in a comparison between RVTRRF and RVT, RVTRRF was found to explain significantly greater variance than RVT in an average of 26% of voxels, whereas RVT explained significantly greater variance than RVTRRF in an average of 4.4% of voxels (p<0.05). This, along with data shown in Table 4, indicates that RVTRRF explained more variance in the resting state than RVT.
A group-level random-effects map of respiratory (PETCO2 and RVTRRF together) contributions to the BOLD signal is shown in Fig. 7. Significant effects were found bilaterally in many regions, consistent wit h (Birn et al., 2006; Wise et al., 2004); here, regions with the strongest effects included cerebellum, insula/orbitofrontal cortex, putamen, caudate, superior temporal gyri, supramarginal gyri, anterior cingulate cortex, mid-cingulate cortex, supplementary motor area, precuneus, and dorsolateral/ventrolateral prefrontal cortex.
As shown in Table 7, the average time delay between PETCO2 and BOLD signal change was 10 ± 1.6 s, consistent with previous reports ((Wise et al., 2007) observed 12 ± 5 s). Time delays between RVTRRF and the BOLD signal were centered near 0 (0.06 ± 2.8 s), which is expected since the derivation of the RRF already accounted for the approximate time delay between RVT and the BOLD signal response. The latency maps for PETCO2 and RVTRRF were significantly correlated (Fig. 8, Table 7).
Regional variability in latency was present within each subject (Fig. 8) as well as across the group (Fig. 9). Figure 9 shows a group-averaged spatial map of the relative latency of BOLD signals to PETCO2; the group map for RVTRRF latency was qualitatively similar and not shown here. While a group-averaged spatial latency map, particularly over a small sample size, can obscure a number of potentially interesting latency relationships, some trends are evident; for example, the cerebellum exhibited a consistently high latency.
While isolated regions with particularly high or low latencies may represent spurious correlations due to noise, we postulate that smoothly-varying latency differences on the order of several seconds may reflect actual differences in vascular response delays. Because a prior study found that a breath-holding (BH) task could be used to query regional differences in the delay of the hemodynamic response uncoupled from CMRO2 changes (Chang et al., 2008), one subject was scanned while performing the BH task described in (Chang et al., 2008) in addition to the resting state scan described above, and latency maps derived from the 2 scans were compared. Here, a time resolution of 0.1 s was used for the cross-correlation analysis. Figure 10 shows a remarkably strong agreement between the BH latency map and the resting-state PETCO2 latency map (r=0.75); the RVTRRF latency map obtained from the resting scan was also consistent with the BH latency map (r=0.75). Furthermore, there is a resemblance between Fig. 9 and the group average latency map from (Chang et al., 2008).
The present results link and extend previous studies which have separately shown that PETCO2 (Wise et al., 2004) and a measurement of respiration volume per time using a standard pneumatic belt (Birn et al., 2006; Birn et al., 2008) are both robust predictors of resting-state BOLD signal fluctuation. Here, by monitoring both physiological signals simultaneously during an fMRI scan, it is demonstrated that PETCO2 and RVT(RRF) have a highly significant linear relationship and, accordingly, account for a remarkably similar pattern of spatial and temporal BOLD signal variation. However, differences between PETCO2 and RVTRRF were present, and it was found that significantly more BOLD signal variance could be explained by a combination of the two signals.
In the present data, neither PETCO2 nor RVTRRF had uniformly greater correlations with the BOLD signal across all subjects as apparent in Figs. 4 and and5;5; some subjects and regions of the brain were more highly correlated with PETCO2, while others were more highly correlated with RVTRRF. Since PETCO2 is believed to be more closely related to mechanisms that directly initiate BOLD signal changes than RVTRRF, it is perhaps surprising that RVTRRF is sometimes a more accurate predictor of the fMRI time series. However, PETCO2 has the potential for being “noisier” than RVTRRF because both the sampling interval and accuracy of PETCO2 are highly dependent on a subject’s behavior: (1) PETCO2 can only be sampled when a subject exhales, whereas RVT yields continuously-sampled data, and (2) since a nasal cannula was used, the PETCO2 measurement could be inaccurate when a subject does not consistently breathe through his/her nose, despite being instructed to do so. On the other hand, cases where RVTRRF yielded a lower fit than PETCO2 may have arisen not because RVT is a less direct measure of relevant breathing processes, but rather because the RRF – being a fixed function derived from previous studies – is less appropriate for certain subjects/brain regions, and a subject- or voxel-specific impulse response may have improved the fit. In addition, the present study used a single respiration belt placed around the upper abdomen; this can be problematic if a subject breathes primarily through the rib cage, or varies the relative contributions of rib cage and diaphragmatic breathing over the scan. The use of 2 measurement belts, one around the abdomen and one around the rib cage, may provide a more accurate estimate of tidal volume (and hence RVT) (Binks et al., 2007; Konno and Mead, 1967).
Across subjects, there was variability in the degree of BOLD signal variance explained by respiratory processes; as can be seen in Fig. 4, some subjects (Subjects 1, 3, and 5) exhibited a much greater effect across the brain for both RVTRRF and PETCO2 than others (Subjects 2, 4, 6, and 7). In an attempt to gain insight into the possible causes of inter-subject variability, several factors were examined: (1) standard deviation of the PETCO2 signal, (2) standard deviation of the PETO2 signal, (3) standard deviation of the RVT signal, (4) motion (drift and RMS), and (5) amplitude of brain signals (average voxel standard deviation across all gray-matter voxels). It was assessed whether any of these factors were correlated with how much overall variance was explained by PETCO2 (or RVTRRF) across the brain. For each of the 4 factors, a correlation coefficient was computed with each of the 3 following metrics of overall variance explained: (1) the number of voxels having Z>5.3, normalized by the total number of voxels in the brain (in order to control for differences in brain size), (2) the mean Z-score across voxels having Z>5.3, and (3) the product of (1) and (2). No significant relationships were found for any of the 4 factors, using any of the 3 metrics of overall variance explained. Therefore, it is possible that individual differences in cardiopulmonary physiology or neurovascular reactivity, as well as the magnitude of neural and/or non-neural noise sources, may play a role. Another possibility is that the dynamics of other physiological processes, such as heart rate and arterial blood pressure, must be considered in order to explain how strongly respiratory effects correlate with BOLD signal time series.
Given that respiratory processes can contribute substantially to the measured BOLD signal time series, removing their effects is an important consideration for fMRI studies of neural function. One must be cautious, however, in studies where physiology is likely to be modulated by a task (e.g. the presentation of emotional or painful stimuli); if physiological noise is task-correlated, disentangling the neural and physiological components can be difficult or impossible. Although the BOLD response to respiration has a different transfer function than that of neural activation (the former reaches a maximum (negative) value at around 15 s, while the latter peaks at approximately 6 s and has a smaller width), breathing-related changes that are out of phase with the task may still appear to be task correlated, particularly in block-design experiments.
CO2 is a powerful modulator of CBF via vasodilation. An increase in blood gas CO2 causes vasodilation and increased CBF, leading to decreased concentrations of deoxyhemoglobin and thus an increase in T2* (BOLD contrast). Conversely, decreases in CO2 cause vasoconstriction and decreased CBF/BOLD signal intensity. Indeed, a simple time-shift of PETCO2 time series was found to predict a significant component of the BOLD signal time series. The PETCO2 fluctuations observed here were largely the immediate effects of respiratory changes, as RVT was negatively correlated with PETCO2 after a time delay of approximately 5 s (1–2 breaths), due to the impact of ventilation on PaCO2. It was further observed that convolving RVT with the RRF impulse response enhanced the cross-correlation with PETCO2. Thus, one way to frame the influence of respiration on the brain is that RVT → CO2 changes → BOLD signal changes, and that the (empirically-derived) RRF essentially captures the transfer function between RVT and CO2, as well the transfer function between the CO2 modulation and the BOLD signal.
It is possible to further describe spontaneous respiration-induced CO2 and BOLD signal changes in terms of a model based on breath holding. It is notable that the RRF (Fig. 1) resembles the early phases of the BOLD response to inspiratory breath holding (Nakada et al., 2001; Thomason et al., 2005), which is marked by (1) an initial positive deflection, followed by (2) a larger negative deflection and (3) a recovery to baseline (and possibly a further increase beyond baseline, if a person continues to hold his breath). This breath-holding response has been suggested to result from the following events (Nakada et al., 2001; Thomason et al., 2005):
Although breath holding typically involves longer and deeper breaths than those occurring naturally in resting state, and although the RRF was first derived by averaging the BOLD response to a series of cued deep breaths (Birn et al., 2008), the fact that a similar response shape maps RVT to the BOLD signal in resting-state data is striking and may indicate the occurrence of similar processes. The variability of heart rate in accordance with the breathing cycle (faster during inspiration, slower during expiration), occurs commonly during spontaneous breathing, and is a phenomenon known as respiratory sinus arrhythmia (Berne and Levy, 1993). Here, in at least one subject, a causal influence from RVT to heart rate was observed (data not shown). While heart rate has been shown to relate to resting-state BOLD fluctuations (Chang et al., 2009; de Munck et al., 2008; Shmueli et al., 2007), the mechanism by which cardiac events translate into BOLD signal change, or modulate respiratory-related BOLD responses, remains to be further elucidated.
Lastly, while we have focused on a forward relationship driven by respiration (Δrespiration → ΔCO2 → ΔBOLD), other autonomic control mechanisms may also operate in the reverse direction. For instance, changes in respiration can result from changes in CO2 via chemoreceptor regulation (Modarreszadeh and Bruce, 1994; Van den Aardweg and Karemaker, 2002), as well as from neural activity due to automatic or voluntary control of breathing (Berne and Levy, 1993; McKay et al., 2003). However, the exploration of such feedback mechanisms is beyond the scope of the current study, as the focus was instead on how changes in breathing pattern causally impact BOLD signal fluctuations as a source of noise. While a number of factors can influence PETCO2 readings, it was observed (Fig. 2) that a large fraction of the resting PETCO2 variability could be attributed to changes in breathing patterns occurring several seconds earlier.
The time delays at which PETCO2 and RVTRRF correlated with voxels in the brain were consistent with one another and revealed patterns of spatial variability. It was furthermore shown that the PETCO2 (and RVTRRF) maps agreed strongly with latency maps derived from a breath-holding task.
Measuring and correcting for vascular delay is an important consideration when performing analyses of functional connectivity/causality as well as event-related activation, as intrinsic vascular delays can mask true directions of neural causality and create spurious or diminished activation (Chang et al., 2008). The correspondence between the latency map from the BH task and those from spontaneous PETCO2 and respiratory-belt derived measures indicates that intrinsic vascular delays may be reliably measured from resting-state respiration data alone, though perhaps with less efficiency and reliability than controlled-breathing tasks. It may also serve as an indication that the physiological processes governing the BOLD response to BH and to spontaneous breathing changes are largely similar.
Measurements of respiration using a standard pneumatic belt, when appropriately transformed (RVTRRF), yielded information about the resting-state BOLD signal that was comparable to that of end-tidal gas monitoring. Thus, to a large degree, pneumatic belt measurements can be utilized as a surrogate for end-tidal gas monitoring. However, the overlap between the two signals was not complete, and the variability across subjects in terms of which of the measurements (PETCO2 or RVTRRF) explained greater BOLD signal variance suggests that, when practical, monitoring both signals is preferable. Finally, the latency at which PETCO2 and respiration-belt measures affected the BOLD signal were consistent and may reveal the spatial distribution of vascular response delays, suggesting that this information may prove useful for calibrating hemodynamic response timing across the brain for activation and connectivity studies when an additional calibration scan cannot be obtained.
This research was supported by NIH grants P41-RR09784 to GHG and F31-AG032168 to CC. We extend special thanks to Sean Mackey for assistance with the capnography instrument, to Jarrett Rosenberg for helpful advice concerning statistics, to all of our participants for volunteering for this study, and to two anonymous reviewers for suggestions that have greatly improved this manuscript.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errorsmaybe discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.