|Home | About | Journals | Submit | Contact Us | Français|
It has previously been shown that low-frequency fluctuations in both respiratory volume and cardiac rate can induce changes in the blood-oxygen level dependent (BOLD) signal. Such physiological noise can obscure the detection of neural activation using fMRI, and it is therefore important to model and remove its effects. While a hemodynamic response function relating respiratory variation (RV) and the BOLD signal has been described (Birn et al., 2008b), no such mapping for heart rate (HR) has been proposed. In the current study, we simultaneously deconvolve the effects of RV and HR from resting state fMRI. We demonstrate that a convolution model including RV and HR can explain significantly more variance in gray matter BOLD signal than a model that includes RV alone, and propose an average HR response function that well characterizes our subject population. We find that the voxel-wise morphology of the deconvolved RV responses is preserved when HR is included in the model, and is adequately modeled by Birn et al.’s previously-described respiration response function. We further demonstrate that modeling out RV and HR can significantly alter functional connectivity maps of the default-mode network.
Functional neuroimaging using MRI (fMRI) relies on the use of blood-oxygen level dependent (BOLD) contrast to depict brain regions that respond to task-induced activation or are functionally connected to other regions (Bandettini et al., 1992; Biswal et al., 1995; Kwong et al., 1992; Ogawa et al., 1992). BOLD contrast results from hemodynamically-driven changes in tissue and vessel oxygenation and is therefore an indirect measure of cerebral metabolism. Unfortunately, physiological processes such as cardiac pulsatility and respiration can also cause changes in cerebral blood flow, thereby inducing substantial fluctuations in the BOLD signal that may confound inferences made about neural processing from analyses of BOLD signals.
Pulsatility of blood and cerebrospinal fluid (CSF) due to cardiovascular processes causes artifacts that tend to spatially localize near ventricles, sulci, and large vessels (Dagli et al., 1999; Glover et al., 2000). Respiration may be accompanied by bulk motion of the head as well as modulation of the magnetic field by thoracic and abdominal movement, and the noise induced in fMRI is more spatially global (Glover et al., 2000). Accordingly, a number of methods have been developed to de-noise fMRI time series by filtering out signals that are time-locked to the cardiac and respiratory phase waveforms, measured by a photoplethysmograph and pneumatic belt, respectively (Deckers et al., 2006; Glover et al., 2000; Hu et al., 1995).
Breathing can also cause a different form of BOLD contrast, thought to result from modulation of blood flow and CO2 in the brain in the presence of ongoing basal metabolism, and corresponding vasomotor regulation (Birn et al., 2006; Corfield et al., 2001; Kastrup et al., 1999a; Kastrup et al., 1999b; Kastrup et al., 1999c; Kastrup et al., 1998; Liu et al., 2002; Nakada et al., 2001). Subtle variations in breathing depth and rate that occur naturally during rest can therefore account for a significant amount of variance in the BOLD signal which, importantly, affects widespread regions of gray matter (Birn et al., 2006; Wise et al., 2004). These low-frequency variations in respiration volume (RV) are especially problematic for studies of task-free resting state, as their spectra overlap with the frequency range of functionally connected networks (<0.1 Hz) (Cordes et al., 2001). Indeed, including RV as a nuisance covariate in a regression model can alter functional connectivity maps of the default-mode network (Birn et al., 2006). Birn et al. further showed that the linear transfer function mapping between RV and the BOLD signal is well modeled by a biphasic curve with a predominantly negative deflection, having an overall duration of approximately 30 sec (Birn et al., 2008b)
A recent study suggested that heart rate (HR) fluctuations may be another source of resting state BOLD signal variance (Shmueli et al., 2007). By including time-shifted HR time series in a general linear model, Shmueli et al. found they explained an aditional 1% of BOLD signal variance beyond RV and RETROICOR regressors. The brain regions in which HR explained additional variance were not concentrated entirely around large vessels, but included gray matter and were sometimes co-localized with regions showing significant correlations with RV. In addition, they observed that HR was negatively correlated with the BOLD signal time series at time lags ranging from 6-12 sec, and positively correlated at time lags of 30-42 sec. This observation indicates the possibility of a more complex temporal relationship between HR and the BOLD signal than is described by cross-correlation. To date, however, a cardiac-related hemodynamic response function has not been studied or even proposed.
In the present study, we employ a linear systems model to relate both RV and HR fluctuations to components of the BOLD signal time series. RV and HR impulse responses are simultaneously deconvolved on a voxel-wise basis using one session of resting state data, and their predictive power is evaluated using a separate session of resting state data from the same subject. One primary aim is to determine whether a convolution model that includes both RV and HR can explain significantly more variance than a model that includes RV alone. Allowing HR to enter the model through a convolution, rather than time-shifted correlations, permits the discovery of a more flexible and accurate mapping between HR and the BOLD signal.
A second aim is to characterize both the RV and HR impulse responses resulting from the simultaneous deconvolution. Even if the inclusion of HR explains significantly more variance, it is not known whether the nature of the mapping varies widely across the affected regions of the brain, or whether a single average response can serve as a representative mapping for most voxels. The deconvolved RV impulse response is also of interest; although an average RV impulse response has been characterized (Birn et al., 2008b), it is not known whether interactions between respiratory and cardiac processes would result in a regionally diverse RV impulse response when HR is also included in the model.
A third aim is to examine the impact of our model’s RV and HR corrections on functional connectivity maps of one particular resting state network, the default-mode network (DMN). The DMN comprises a set of regions that exhibit low-frequency correlated signals in task-free resting state (Greicius and Menon, 2004; Raichle et al., 2001), and collectively down-regulate during a wide range of cognitively demanding tasks (Binder et al., 1999; Gusnard et al., 2001; Mazoyer et al., 2001; McKiernan et al., 2003; Shulman et al., 1997). Quantitative measurement of connectivity in networks such as the DMN is increasingly employed to draw inferences about behavior (Clare Kelly et al., 2008; Daselaar et al., 2004; Hampson et al., 2006), development (Damoiseaux et al., 2007; Thomason et al., 2008), and dysfunction (Garrity et al., 2007; Greicius et al., 2007; Greicius et al., 2004; Uddin et al., 2008); therefore, the reduction of physiological noise that might confound delineation of the network’s neural characteristics is critical.
Participants included 10 right-handed, healthy adults, including 3 females (mean age = 31.4 years, SD=13.4). All subjects provided written, informed consent, and all protocols were approved by the Stanford Institutional Review Board.
Subjects underwent 2 scans during which no intentional task was performed (“Rest1” and “Rest2”), with respective durations of 12 min and 8 min. Subjects were instructed to relax and close their eyes while remaining awake. Between the 2 resting state scans, subjects performed a 10 min event-related Sternberg working memory task, consisting of 3 back-to-back 5 sec trials (0.5 sec encoding, 3 sec maintenance, 1.5 s probe) followed by 45 sec of fixation, i.e. having a mean ISI of 60 sec. Encoding stimuli consisted of 4 uppercase letters in a cross-like configuration around the center of the screen, and the probe stimulus was a single lowercase letter presented at the center of the screen. This task was used only to localize subject-specific seed regions of interest (ROI) for the resting state connectivity analysis (described below). No further analysis was performed on these data.
Magnetic resonance imaging was performed on a 3.0-T whole-body scanner (Signa, rev 12M5, GE Healthcare Systems, Milwaukee, WI) using a custom quadrature birdcage 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, 1-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 misregistration occurs from off-resonance in spiral imaging (Pfeuffer et al., 2002).
Heart rate and respiration were monitored at 40 samples/sec using the scanner’s built in photoplethysmograph placed on the left hand index finger and a pneumatic respiratory belt strapped around the upper abdomen, respectively. A file containing cardiac trigger times and respiratory waveforms was generated for each scan by the scanner’s software. Values in the respiratory waveform were converted to a percentage of the full scale (difference between the maximum and minimum belt positions measured over the scan). Only the fractional variations in the waveform, rather than the absolute amplitude values, are of importance in the current study.
Motion parameters were calculated using methods described in (Friston et al., 1996). Coregistration, however, was not performed since (1) motion was expected to be minimal because a bite-bar was used, (2) coregistration causes unintended smoothing across voxels, which would interfere with our voxel-wise analysis, and (3) estimation of coregistration parameters can be biased by activation in tasks that evoke large, widespread signal changes (Freire and Mangin, 2001). Therefore, it was important to verify that intra- and inter-scan motion was minimal. We calculated the maximum peak-to-peak excursion, root mean square (RMS) fluctuation, and task-correlated motion for the 3 translational and 3 rotational motion parameter time series within each scan (rotations were converted to worst-case translations by multiplying by 65mm, an average head radius (Thomason and Glover, 2008)). Summary statistics are reported as the maximum of these values over the 6 axes of motion.
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=5mm), and removal of linear and quadratic temporal trends. Spatial normalization to a standard template was not performed, to avoid blurring; all computations were done in the original subject-space, and voxels were maintained in their original dimensions (3.4375mm × 3.4375mm × 4mm). The first 7 temporal frames were discarded to allow the MR signal to equilibrate.
Following pre-processing, data were corrected for cardiac pulsatility and respiratory motion artifacts using RETROICOR (Glover et al., 2000). Thus, the RV and HR results presented in the current study represent noise contributions beyond those merely synchronous with the cardiac and respiratory cycles.
The respiratory waveform recorded by the pneumatic belt is related in a complex manner to the subject’s tidal breathing volume and other pulmonary characteristics. Nevertheless, a measure that is loosely associated with tidal volume, and hypothetically BOLD signal modulation, can be extracted from this waveform in various ways. Birn et al. calculated the respiration volume per unit time (RVT), defined as the difference between the maximum and minimum belt positions at the peaks of inspiration and expiration, divided by the time between the peaks of inspiration (Birn et al., 2006). Here, we propose a measurement (referred to as RV) that is based on computing the standard deviation of the respiratory waveform on a sliding window of 3 TRs (6 sec) centered at each desired TR sampling point. Specifically, the value of the RV time series assigned to the kth TR was computed by taking the standard deviation of the raw respiratory waveform over the 6 sec time interval defined by the (k-1)th, kth, and (k+1)th TRs. Thus, RV(k) is essentially a sliding-window measure related to the inspired volume over time. This measure is simpler and more robust to noise than RVT, as it calculates the RMS average fluctuation over a window rather than taking a single peak-to-valley difference, and does not rely on the accuracy of peak detection required for breath-to-breath computations. However, because differences between RV and RVT are expected, we performed comparisons between the two waveforms (described below). The HR(k) time series was computed by averaging the time differences between pairs of adjacent ECG triggers contained in the 6 sec window defined by the (k-1)th, kth, and (k+1)th TRs, and dividing the result into 60 to convert it to units of beats-per-minute.
The time series of a voxel (y) is modeled as the sum of RV convolved with an unknown RV filter (hr) and HR convolved with an unknown HR filter (hh), plus a noise term (ε):
where hr ~ N(0,Kr), hh ~ N(0,Kh), and .Xr and Xh refer to the convolution matrices defined by RV and HR, respectively, and Kr, Kh and will be defined below.
By defining the filters hr and hh to be Gaussian processes (which here, since time is discretized, can be considered Gaussian random vectors), we enforce temporal smoothness while allowing greater flexibility in shape compared to approaches such as constraining the filters to lie in the span of pre-specified basis functions; an appropriate basis set for describing the RV and HR impulse responses was not known a priori. The use of Gaussian process priors for deconvolution has previously been applied in the fMRI literature in the context of estimating the hemodynamic response function (HRF) to task activation (Casanova et al., 2008; Goutte et al., 2000; Marrelec et al., 2003a; Marrelec et al., 2003b, 2004), where it was demonstrated to better capture features of the HRF, such as the post-stimulus undershoot, than sets of Gamma and Gaussian basis functions (Goutte et al., 2000).
The model (1), which will be referred to as the RVHR model, can be written more compactly as y = Xh +ε, where X =[Xr,Xh ] and h = [hr,hh ]T. Then h ~ N(0,K) with , and the maximum a posteriori (MAP) estimate of h is
(see Appendix A). The covariance matrices Kr and Kh describe the degree of correlation between points in hr and hh as a function of their distance. We let
This form, known as the squared exponential kernel, is a standard choice in the use of Gaussian processes for regression (Rasmussen and Williams, 2006). The length scale l governs the degree of smoothness imposed on the deconvolved filter (increasing l will produce more slowly-varying filters) and the kernel variance regulates the distance from which values of h depart from its mean (which is 0 in this case).
One might choose to optimize all 3 hyperparameters (l, , and ) at each voxel by maximizing the associated likelihood function. However, to reduce the degrees of freedom and potential for overfitting (and to increase computational efficiency), we fixed l = 2 and for all voxels; at the ith voxel was then set equal to the sample variance of the ith voxel time series. These values were selected from preliminary experiments in which optimization of the likelihood function over all 3 hyperparamters was performed using a nonlinear conjugate gradient method (using a MATLAB implementation by (Rasmussen, 2006)) where it was observed that the distributions of l and estimates were centered at these values. When fixing l and the estimate was proportional (and close) to the sample variance of y. We found that the goodness of fit and resulting deconvolved filters were similar when fixing the hyperparameters or when allowing all 3 to vary, although fixing the hyperparameters (particularly l) tended to produce more slowly-varying filters.
Both hr and hh were additionally constrained to be 0 at their beginning and ending, as we assume that the hemodynamic response to sudden fluctuations in either RV and HR has a finite rise time and that the influence ultimately decays away (see Appendix B for details of the implementation). Both hr and hh were assigned durations of 30 sec (15 points).
A reduced model in which the HR input is absent was also implemented. In this RV model,
where the assumptions and estimation of r are as described in the RVHR model. We also consider a model in which the RV filter at every voxel was defined to be the “respiratory response function” (RRF) defined by (Birn et al., 2008b); we refer to this as the RRF model.
For each subject, voxel-wise deconvolution of hr, hhand r was performed on the Rest1 scan; the generalizability of the deconvolved filters, as well as model comparison, was evaluated by applying the models to the Rest2 scan. For the RVHR model, RV and HR time series were extracted from the Rest2 physiological data and convolved with the RV and HR voxel-wise impulse responses obtained from Rest1, forming 2 unique covariates for each voxel. A linear regression against these covariates was performed for each voxel’s Rest2 time series, providing both a fitted response and a residual error value as well as a residual time series, which we define to be the “corrected” signal. The same was performed for the RV model, in which the voxel-wise impulse responses from the RV model obtained from Rest1 were convolved with the RV waveform from Rest2 to obtain a single unique covariate for each voxel. The same was also performed using the RRF as the impulse response for every voxel (RRF model). The significance of variance explained by each model, as well as by the RV and HR components of the RVHR model, was computed using an F-test.
The variance explained by the RVHR, RV, and RRF models were compared for 3 pairs of models: (1) RVHR>RV, (2) RVHR>RRF, and (3) RV>RRF. The aim of the first test was to determine whether modeling RV and HR has a significant effect beyond modeling RV alone; the aim of the second test was to compare the full model to the current standard (modeling only RV, using a canonical response function); the aim of the third test was to determine whether using voxel-wise RV impulse response functions is more effective than using a canonical RV impulse response. Comparisons (1) and (2) were computed using an F-test; since the two models in comparison (3) are of the same complexity, the difference between correlation coefficients from the 2 models was tested for significance. This was performed by transforming the r values to Fisher Z values, normalizing by (here with n = 233 frames), and assessing the significance under the standard normal distribution (Dowdy and Wearden, 1991).
Functional connectivity maps of the DMN were compared before and after correction with the RRF, RV, and RVHR models. A DMN functional connectivity map was obtained for individual subjects by extracting the average time series from a seed ROI in the posterior cingulate cortex (PCC)/precuneus, a central node of the DMN (Greicius et al., 2003; Greicius and Menon, 2004; Raichle et al., 2001) and correlating it against every voxel in the brain. The seed ROI was defined as the intersection of a template anatomical PCC/precuneus (MarsBar AAL atlas; http://www.sourceforge.org/marsbar) with the subject’s deactivation map in the WM task (thresholded at r<-0.2). The template PCC/precuneus ROI was normalized to the subject’s average functional image over the Rest1 scan using SPM5 (http://www.fil.ion.ucl.ac.uk/spm). The procedure of extracting a seed time series and correlating it with all voxels was performed on the Rest2 datasets before any physiological noise correction was applied, and repeated after each of the above corrections (RV, HR, RRF) were applied, with and without RETROICOR. Group maps of the DMN were computed by normalizing Fisher-Z transformed correlation maps to a standard EPI template and entering them into a random-effects analysis using SPM5 (http://www.fil.ion.ucl.ac.uk/spm).
RV and RVT are computed in slightly different ways from the respiration belt measurements. While RV represents the RMS average inspired volume over a 6-sec sliding window, RVT considers a shorter interval (breath-to-breath, which tends to be 3-4 sec) and accounts more explicitly for variations in breathing rate by normalizing the depth by the breath-to-breath time interval. While we expected that differences between the RV and RVT waveforms would be minor, it was important to examine the degree to which results depend on the choice of one or the other. Therefore, the variance explained by each model (as described in the above section, “Evaluation”) was also computed using RVT in place of RV, and the resulting average deconvolved filters were compared.
To gauge the influence of RETROICOR as a pre-processing step, the RVHR model was also fit to the data without first performing RETROICOR. DMN connectivity maps were then compared for 3 cases: correction with RETROICOR only, RVHR only, and both RETROICOR and RVHR. Changes were quantified by counting the number of voxels having significant (p<0.05) correlation with the PCC seed ROI for each correction, and converting the count into a percentage (with respect to the number of significant voxels in the uncorrected – that is, without RETROICOR or RVHR – connectivity maps).
The estimated motion parameters are summarized in Table 1. Motion was minimal, as subjects exhibited a mean drift of 1.15 mm (<1 voxel) across the 2 resting state sessions combined.
Summary statistics for the RV and HR measures are shown for each subject in Table 2A. Over all subjects and both resting state scans, the mean RV fluctuation was 16.6 ± 4.4%, while HR fluctuated about 61.2 ± 3.1 beats per minute. RV and HR were only mildly correlated (Table 2B).
Table 3 indicates the fraction of voxels (relative to the whole brain) having significant (p<0.0001) variance explained by the RVHR, RV, and RRF models in the Rest2 scan. The percentage of signal variance explained, averaged over the set of significant voxels for each model, is also provided. For all but one subject, the RVHR model explained significant variance over a greater spatial extent than the RV and RRF models, and accounted for a larger mean percentage of the signal variance.
Maps depicting the percent signal variance explained at each voxel for the RRF, RV, and RVHR models in the Rest2 scan are shown for 3 subjects in Figure 2; the RVHR model is further decomposed into separate RV and HR components in Figure 3. Regions correlated with HR tended to comprise gray matter, and were often disjoint from regions correlated with RV.
For all subjects, the RVHR model explained significant additional variance over the RV or RRF models over some extent of the brain. Percentages of subjects’ brains in which the model comparisons (1) RVHR>RV, (2) RVHR>RRF, and (3) RV>RRF were significant are presented in Table 4.
Correcting for RV and HR tended to decrease the overall spatial extent of DMN connectivity. Figure 4 shows the functional connectivity of the PCC for 2 subjects without correction (A), and with correction using the RV model (B) and RVHR model (C); correction using the RRF model was similar to (B). Functional connectivity with the PCC tended to become increasingly focal with correction, though the amount of change due to each type of correction varied across subjects. Significant (p<0.05) reductions in connectivity among the set of initially-connected voxels (i.e. those having r>0.11 (p<0.05) without correction) occurred for an average of 1.0% (SE=0.3%) of voxels after the RRF model correction, 5.6% (SE=4.3%) after the RV model correction, and 7.6% (SE=3.8%) after the RVHR model correction. Respective increases were all below 2%. A comparison of group-level maps before any physiological noise correction (Figure 5A) and after both RETROICOR and RVHR corrections (Figure 5B) show a reduction of apparent false-positives, as well as increased connectivity in areas believed to be canonical nodes of the DMN, such as the medial prefrontal and anterior cingulate cortices.
Average deconvolved HR and RV filters from the RVHR model (hr and hh) were computed for each subject from voxels in which the HR and RV components reached significance (p<0.001) by a t-test on their respective regression coefficients. The grand average was then taken, yielding the filters plotted in Figure 6(A,B). The average RV filter bears a strong resemblance to the RRF (Figure 6C), while the average HR filter exhibits a peak at 4 sec and a dip at 12 seconds. The standard error bars for the average RV and HR filters are comparable and small, indicating consistency across subjects. To further query the spatial properties of the HR filters, k-means clustering was performed on the voxel-wise HR filters pooled over all subjects, using 6 clusters and a similarity metric of correlation. Figure 7 shows the mean HR filter from each cluster, and the spatial locations of the clusters for 2 subjects. The clusters encompassing gray matter had mean filters resembling the average HR filter, indicating minor spatial variability.
To examine whether the average HR filter obtained from our subject population is generalizable, we applied it to a set of 3 subjects from a separate study who had each undergone one 8 min resting state scan. An analytic function (CRF(t), for “cardiac response function”) was fit to the curve as the sum of a Gamma and a Gaussian:
(Figure 6D), and this function was used as the HR filter for all voxels. Two models were then evaluated on each of the 3 subjects’ resting state scans: (1) the RRF model, as described previously; and (2) the “RRF-CRF” model, in which the HR time series, convolved with CRF(t), was included as a covariate in the linear model. The RRF was used as the RV filter for each voxel since a separate scan was not available for deconvolving voxel-wise RV responses without bias, and because the RRF has been established both in the present study and in (Birn et al., 2008b) to well model the transfer function between RV and the BOLD signal. Maps of the variance explained by each model are shown for 2 subjects in Figure 8. Including HR convolved with CRF(t) contributed a significant amont of variance in all subjects’ resting state scans (Table 5), suggesting that CRF(t) is perhaps a generalizable transfer function between HR and the BOLD signal.
Strong correlations were obtained between estimates of the RV filter resulting from the RVHR and RV models; that is, voxel-wise estimates of hr and r were highly consistent despite the inclusion of HR in the RVHR model. The mean correlation across the set of voxels in which significant variance (p<0.001) was accounted for by both the RV and RVHR model was r=0.97 (SD=0.05); even when restricting the analysis to the set of voxels havng a significant (p<0.001) HR component in the RVHR model, strong correlations were maintained (r=0.93, SD=0.10). Similar generalization performance between the RV and RRF models reflects the fact that hr and r were also well correlated with the RRF (r=0.74, SD=0.11 between r and the RRF).
The RV and RVT waveforms were correlated (mean r=0.61, SD=0.1 over the Rest2 scans of all 10 subjects). Minor differences were found in the shapes of the average filters when RVT was used instead of RV (Figure 9), though in comparing Figures 9A and and6A,6A, we note that the error bars for each time point were slightly smaller when RV was used. Differences in the percentage variance explained when RVT was used instead of RV are shown in Table 6. While differences were again minor, there was a trend towared greater average variance explained when RV was used.
Figure 10 reveals how the extent of DMN connectivity changed after correction using RETROICOR, RVHR, and both RETROICOR and RVHR. It is apparent that (1) applying RETROICOR alone influenced the outcome of DMN connectivity slightly, (2) connectivity changes due to RVHR were much greater than those induced by applying RETROICOR (p<0.05), and (3) the extent (volume) of change due to applying RETROICOR before RVHR was similar to that of applying RVHR alone (p>0.05; n.s.). Figure 11 illustrates the connectivity maps for one subject.
The present study demonstrates that a linear systems model having both RV and HR inputs can account for substantial fluctuations in the resting state BOLD signal. The RVHR model explained a significantly greater fraction of signal variance beyond the RV model, and over a spatial extent encompassing widespread regions of gray matter. A HR hemodynamic response function (CRF(t)) is proposed, and is shown to adequately characterize the mapping between HR and the BOLD signal for our subject population. Furthermore, results verify that the RV response function (RRF) introduced by Birn et al. represents well the mapping between RV and the BOLD signal, even when simultaneously accounting for HR. Moreover, removing HR and RV components from the BOLD signal using our RVHR model can induce significant changes (region-specific reductions as well as increases) in resting state network connectivity.
The RV and HR filters were deconvolved using Gaussian process priors. This yields smooth (and thus physiologically plausible) response functions, and so presents a suitable approach to the problem of deconvolving physiological contributions to the BOLD signal. This method avoids the noise amplification problem of unconstrained deconvolution which has led others (e.g. (Birn et al., 2008b)) to use less straightforward approaches. In addition, the framework for simultaneous deconvolution employed in the current study can be directly extended to incorporate additional physiological or neural inputs that may be of future interest.
In agreement with Shmueli et al., HR was found to contribute significant additional variance beyond RV in gray matter as well as in regions expected to be more susceptible to cardiac pulsatility, such as CSF and major vessels (Shmueli et al., 2007). Importantly, as illustrated in Figure 3, regions in which the HR component had a significant effect were often disjoint from regions in which the effects of RV (and RRF) were significant, indicating differential contributions of HR and RV effects over space. Variability was observed across subjects with respect to both magnitude and spatial location of significant HR effects, and a group analysis did not reveal any particular anatomical regions that were highly significant for HR. While the mechanism through which HR could modulate the BOLD signal is not well-understood, gray matter correlations with HR could be due in part to neuronal activity linked with changes in levels of arousal.
Deconvolution of the HR response revealed that the majority of responsive voxels exhibited a peak at around 4 sec and an undershoot of approximately equal magnitude at around 12 sec. These features were remarkably consistent across subjects and anatomical regions, and generalized well to a separate population of 3 subjects. The use of this proposed average HR filter (CRF), in conjunction with the average RV filter (RRF) proposed by Birn et al., may thus present an alternative to voxel-wise deconvolution for reducing the influence of HR and RV in the BOLD signal, obviating the need for an additional scan.
However, if obtaining voxel-wise RV and HR filters is desired, appending a 12 min scan in order to implement the RVHR model on other functional scans would be expensive and time-consuming; thus, it is of interest to know whether a shorter scan will suffice. We evaluated the stability of deconvolved HR responses for one subject using only the first 4 min and 8 min of the 12 min Rest1 scan, for voxels in which the RVHR model was significant (p<0.001). As shown in Figure 12, 8 min of data sufficed to yield responses that were highly similar to those from 12 min.
Voxel-wise estimates of the RV impulse response were only slightly more effective (and for 3 subjects, less effective) in accounting for variance in the second resting state scan compared to using the RRF for all voxels. In addition, the amount of spatial ovelap was substantial, consistent with the fact that our deconvolved RV impulse responses were similar in character to the RRF – even when simultaneously accounting for HR, and despite the fact that RV (rather than RVT) was used as the respiration input function. This finding suggests that the RRF is indeed an effective a priori RV response function, which is further demonstrated by our result showing that the RRF correction impacts maps of the DMN, a finding that extends Birn’s results from correction based on cross-correlation (Birn et al., 2006) to use of the RRF (Birn et al., 2008b). Furthermore, convolving RV with the RRF was shown to provide a better fit to resting state data than a simple cross-correlation; this is somewhat different from the results of Birn et al., 2008b, who had found that the RRF accurately modeled BOLD signal changes to cued-breathing paradigms but did not fit resting state data more effectively than cross-correlation
We chose to deconvolve RV and HR simultaneously, rather than HR alone, because (a) RV is known to account for significant portion of the BOLD signal and it is therefore critical to model it, and (b) possible interactions between RV and HR may result in different filters than if either were estimated alone. Respiratory and cardiac processes are intimately related, and both govern cerebral blood flow and oxygenation and hence BOLD signal changes. Yet, we observed strong correlations between RV filters deconvolved using the RVHR and RV models, even for the set of voxels demonstrating a significant HR component. That modeling HR had a negligible impact on the deconvolved RV filters suggests a degree of independence between the way in which RV and HR may propagate into the BOLD signal.
The RVHR corrections altered DMN maps beyond corrections using RV alone, though both models acted primarily to diminish connectivity with the PCC. It is not possible to quantify the goodness of network changes, as there is no ground truth for which regions ought to comprise the network, and it is not clear to what extent previous characterizations of the DMN reflect physiological, rather than neural, connectivity. However, “canonical” nodes of the DMN are believed to be those showing task-correlated deactivation (Greicius et al., 2003; Raichle et al., 2001) and reduction of connectivity outside of these regions is viewed as desirable. In our study, brain regions were not affected uniformly by RV or RVHR correction; with reference to Figure 5, connectivity of canonical nodes such as the medial prefrontal cortex and anterior cingulate cortex was preserved or increased, whereas regions outside of the network tended to decrease. Hence, modeling both RV and HR can potentially reduce the number of false positives in functional connectivity analysis.
Importantly, Figures 4 and and55 also reveal that the spatial extent of RV and HR fluctuations overlaps substantially with at least one of the resting state networks, the DMN. This suggests that data-driven approaches such as ICA will not be able to fully distinguish between neural sources of connectivity and those due to physiological processes such as RV and HR (Birn et al., 2008a), and therefore, that preprocessing resting state data with the RVHR correction would be useful, especially when quantification of network extent is desired.
It is fairly common in the literature to apply a “global signal” correction, in which a whole-brain average signal is regressed out of each voxel time series prior to analyzing functional connectivity (e.g. (Di Martino et al., 2008; Fair et al., 2008; Fox et al., 2005)). While having the desired effect of making networks more focal, the global signal represents a mixture of physiological, neural, and movement-related sources, and regressing it out will unfortunately remove common signal as well as common noise. Large networks whose constituent regions exhibit large BOLD signal changes, such as the DMN, may contribute substantially to the global average. On the other hand, approaches that model and correct for physiological noise based on external measures, such as the respiration and cardiac monitoring, can reduce potential noise in an unbiased way.
Because we aim to model and remove physiological noise due to cardiac and respiratory processes beyond those simply time-locked to the cardiac and respiratory cycles, correction was initially intended to be applied post-RETROICOR. We observed that quantitative changes in the spatial extent of variance explained (and DMN connectivity) when applying the RVHR model with and without RETROICOR were small; however, such a quantification does not characterize any possible spatial dispartity or ’goodness’ of the change, nor does it reveal any overlap between the effects of both corrections. It is interesting that although the signal components treated by RETROICOR and RVHR are in a sense distinct (though having a common physiological origin), the effects from correction using the 2 methods can be similar, as seen in the DMN for the subject in Figure 11. RETROICOR and RVHR produced comparable changes in the DMN. In this subject, as well as several others, applying RETROICOR with RVHR produced greater changes in DMN connectivity (again, mainly reductions) than applying either in isolation.
Our measure of RV consisted of the sliding-window standard deviation of the respiratory belt waveform, and thus relates primarily to changes in inspiration depth over time. While simpler to compute than RVT (no peak detection in the respiration waveform is required), the use of a fixed-length sliding window rather than breath-to-breath intervals is more difficult to interpret due to its complex dependence on the depth, frequency, and phase of respiratory cycle within each window, especially when the breathing rate is slow relative to the window length. However, our RV measure tended to correlate with RVT) in our set of subjects, and the average RV transfer function agreed with the RVT transfer function (RRF) of (Birn et al., 2008b). In addition, a re-analysis of our data using RVT instead of RV yielded minor differences in the results. While it is encouraging that two slightly different quantifications of respiratory variation from a pneumatic belt waveform can explain similar variance in the BOLD signal, the optimal way to relate measurements of chest expansion during breathing to parameters that more directly mediate BOLD signal change, such as end-tidal CO2, requires further investigation.
In summary, a canonical HR response function is proposed that, in conjunction with the RRF and with or without RETROICOR, can result in significant and widespread reductions of variance in resting BOLD contrast timeseries when applied in a linear model. The relative contribution of HR and RV to BOLD fluctuations was not uniform across subjects; sometimes HR was more prominent than RV, and in other cases less so. In general, applying RVHR corrections resulted in diminished spatial extent of the resting state DMN to a greater degree than with RV corrections alone. This suggests that caution should be exercised when employing quantitative analyses of resting state networks, without accounting for physiological noise sources from both cardiac and respiratory functions.
This research was supported by NIH grants P41-RR09784 to GHG, T32-GM063495 to CC, and a Stanford Graduate Fellowship to JPC. We thank all of our participants for volunteering to be scanned for this study, and two anonymous reviewers who provided important suggestions for improving this manuscript.
As described in the text, the model (Eq. 1) can be written compactly as y = Xh+ε, and a Gaussian process prior is placed on h: h ~ N(0,K). In the Bayesian framework, we would like to find h to maximize p(h | y), which is known as the maximum a posteriori (MAP) solution. Following Bayes’ rule, the posterior distribution of over h is given by
where (I is the identity matrix) and the marginal likelihood p(y) is independent of h. Including only the terms that are dependent on h (i.e. the numerator of Eq. 6), we obtain
Since maximing p(h | y) is equivalent to minimizing -log p(h | y), the problem reduces to minimizing
Setting the gradient with respect to h of Eq. 8 to 0 and solving for h yields
which is the result provided in Eq. 2.
When estimating hr and hh, we wish to enforce the sensible constraint that both filters begin and end at 0. Optimization theory offers a method to incorporate this constraint into the MAP estimate of Eq. 9, which is based on the KKT conditions and described in many optimization texts (see, e.g., (Boyd and Vandenberghe, 2004)). A brief overview of our treatment follows.
To implement the constraint that hr and hh (which are of length nr and nh, respectively) each begin and end at 0, we solve the quadratic optimization problem
where S is a 4x(nr+nh) selector matrix consisting of ones at the (1,1), (2, nr), (3, nr +1) and (4, nr +nh) entries and zeros elsewhere, and R = K-1. The objective function in Eq. (10) results from wishing to maximize p (h | y) (Eq. 8). The KKT conditions specify that the optimal solution satisfies
where h* and υ* are the optimal primal and dual variables. The condition (12) results from the fact that the gradient of the Lagrangian must vanish at h* and υ*. Equations (11) and (12) define the linear system
which can be solved efficiently for both h* and υ*.
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 errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.