PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Neuroimage. Author manuscript; available in PMC 2008 August 1.
Published in final edited form as:
PMCID: PMC2214855
NIHMSID: NIHMS27952

A Component Based Noise Correction Method (CompCor) for BOLD and Perfusion Based fMRI

Abstract

A component based method (CompCor) for the reduction of noise in both blood oxygenation level dependent (BOLD) and perfusion-based functional magnetic resonance imaging (fMRI) data is presented. In the proposed method, significant principal components are derived from noise regions-of-interest (ROI) in which the time series data are unlikely to be modulated by neural activity. These components are then included as nuisance parameters within general linear models for BOLD and perfusion-based fMRI time-series data. Two approaches for the determination of the noise ROI are considered. The first method uses high-resolution anatomical data to define a region of interest composed primarily of white matter and cerebrospinal fluid, while the second method defines a region based upon the temporal standard deviation of the time series data. With the application of CompCor, the temporal standard deviation of resting state perfusion and BOLD data in gray matter regions was significantly reduced as compared to either no correction or the application of a previously described retrospective image based correction scheme (RETROICOR). For both functional perfusion and BOLD data, the application of CompCor significantly increased the number of activated voxels as compared to no correction. In addition, for functional BOLD data, there were significantly more activated voxels detected with CompCor as compared to RETROICOR. In comparison to RETROICOR, CompCor has the advantage of not requiring external monitoring of physiological fluctuations.

Introduction

Over the last decade, blood oxygenation level dependent (BOLD) and perfusion-based functional magnetic resonance imaging (fMRI) have become indispensable tools for studies of the working brain. When utilized together, the BOLD and perfusion signals can provide a quantitative understanding of the metabolic response to neural activity and provide insight into neurovascular coupling mechanisms (Hoge et al. 1999). However, as the fMRI community has moved to higher field strengths, physiological noise has become an increasingly important confound limiting the sensitivity and the application of fMRI studies (Kruger and Glover 2001; Liu et al. 2006).

Physiological fluctuations have been shown to be a significant source of noise in BOLD fMRI experiments, with even a greater effect in perfusion-based fMRI utilizing arterial spin labeling (ASL) techniques (Kruger and Glover 2001; Restom et al. 2006). Physiological sources of noise primarily include cardiac pulsations and respiratory-induced modulations of the main magnetic field. Additional sources include blood flow changes coupled to end-tidal C02 and vasomotion occurring at 0.1 Hz (Hu et al. 1995; Dagli et al. 1999; Glover et al. 2000a).

Approaches to removing cardiac and respiratory related-noise include temporal filtering, image based retrospective correction (RETROICOR), k-space based correction (RETROKCOR) and navigator echo based correction (DORK) (Hu et al. 1995; Biswal et al. 1996; Josephs et al. 1997; Glover et al. 2000a; Pfeuffer et al. 2002). More recently, RETROICOR has been extended to a general linear model (GLM) framework (Lund et al. 2006) and modified for use in ASL studies (Restom et al. 2006). A recent adaptation for BOLD based imaging employs additional regressors describing variations in respiratory volume (Birn et al. 2006).

An alternate approach to the use of external measures of physiological activity or specially modified pulse sequences is to globally subtract average time-courses from regions unlikely to be associated with neural activity, such as ventricles and large vessels (Petersen et al. 1998; Lund and Hanson 2001). However, since this technique removes only the average time-series, it is unable to account for voxel-specific phase differences in the noise due to physiological fluctuations. Additionally, component based techniques, utilizing independent component analysis (ICA) or principal components analysis (PCA), have shown potential in identifying spatial and temporal patterns of structured noise (Thomas et al. 2002; McKeown et al. 2003; Beckmann and Smith 2004). However, the utility of component based methods has been limited to BOLD studies with sampling times short enough to clearly differentiate cardiac and respiratory elements from evoked responses (Thomas et al. 2002), in which case a temporal band pass filter would be adequate for noise removal.

In this paper we present and characterize a novel component based method (CompCor) for the correction of physiological noise in BOLD and perfusion-based fMRI. We show that principal components derived from noise regions-of-interest (ROI) are able to accurately describe physiological noise processes in gray matter regions. In our presentation we investigate the use of two different methods for determining noise ROIs. The first method uses anatomical data to identify white matter and CSF voxels, while the second method uses the temporal standard deviation (tSTD) of the time-series data to identify voxels dominated by physiological noise. We show that the use of principal components derived from a noise ROI as nuisance regressors in a GLM of the fMRI signal can significantly reduce the temporal standard deviation in resting state scans and increase the sensitivity of functional BOLD and perfusion-based studies.

Theory

CompCor Algorithm

The underlying assumption in the CompCor algorithm is that signal from a noise ROI can be used to accurately model physiological fluctuations in gray matter regions. The term “noise ROI” refers to areas (e.g. white matter, ventricles, large vessels) in which temporal fluctuations are unlikely to be modulated by neural activity and are primarily a reflection of physiological noise. The ability to model gray matter physiological noise elements is then predicated on the similarity between physiological fluctuations in the noise ROI and gray matter. A principal components analysis (PCA) is used to compactly characterize the time-series data from the noise ROI. Significant principal components are then introduced as covariates in a general linear model (GLM) as an estimate of the physiological noise signal space.

In this paper we investigate the use of two methods for determining the noise ROI. The first method uses anatomical data to identify voxels that consist primarily of either white matter or cerebrospinal fluid (CSF). Since neural activation is localized to gray matter, fluctuations in white matter and CSF regions should primarily reflect signals of non-neural origin, such as cardiac and respiratory fluctuations.

In the second method, voxels with high temporal standard deviation (tSTD) are used to define a noise ROI. This approach is based on previous preliminary work by Lund et al. (2001) in which areas of high temporal standard deviation were found to correspond to ventricles, edge-regions, and vessels. The advantage of this method is that it utilizes the time-series data to identify a noise ROI without the need for an anatomical scan.

General Linear Model for ASL and BOLD

The general linear model for the BOLD-weighted data can be written as

equation M1
[1]

where b represents the measured BOLD data, Xh represents the stimulus response where X is a N×k design matrix and h is a k×1 vector of hemodynamic parameters. In the case of a block design, X reduces to a vector containing the smoothed stimulus pattern and h reduces to a scalar representing the unknown amplitude. Nuisance parameters are integrated in Sd, where S is a N×l matrix comprised of l nuisance model functions and d is a l×1 vector of nuisance parameters. We have also added physiological noise terms Pc where P is a N×m matrix containing m regressors and c represents the unknown regressor weights. Finally, n represents the additive noise term.

A general linear model (GLM) for ASL data in gray matter can be written as

equation M2
[2]

where p is the acquired raw data representing interleaved tag and control images. In this model, the term modeling perfusion Xhperf is modulated by a diagonal matrix, M, consisting of alternating −1's and 1's for the tag and control images, respectively (Mumford et al. 2006; Restom et al. 2006). The term XhBOLD models a BOLD weighted static tissue component.

In a noise ROI, where we expect no stimulus-related response, the GLM reduces to b = Sd + Pc + n or p = Sd + Pc + n for BOLD and perfusion data, respectively. An assumption of CompCor is that the space spanned by the significant principal components from the noise ROI can be used to estimate the space spanned by the columns of the physiological regressor matrix P.

Figure 1 depicts the basic algorithm in which significant principal components from a noise ROI are used to form a physiological noise matrix Pest which is then used as an estimate of P in the GLM. For functional studies, an added processing step is included to minimize the possibility of including stimulus related fluctuations. In this step, a preliminary GLM analysis, using the appropriate design matrix X , is used to exclude voxels from the noise ROI with a calculated p-value less than 0.2. Additional details are provided in the Methods section.

Figure 1
Schematic of the CompCor algorithm in which significant principal components derived from time-series data within noise regions-of-interest (nROI) are used to form an estimate Pest of the physiological noise matrix P. Incorporation of Pest into the general ...

Methods

Experimental Protocol

Ten healthy subjects (ages 23 to 39) participated in the study after giving informed consent. Each subject viewed one periodic single trial visual stimulus consisting of a 20-second initial “off” period followed by 5 cycles of a 4 second “on” period and a 40 second “off” period. In addition to a periodic design, each subject viewed one block design consisting of 4 cycles of a 20 second “on” period and a 40 second “off” period. During the “on” periods, a full-field, full contrast radial 8 Hz flickering checkerboard was displayed, while the “off” periods consisted of a gray background of luminance equal to the average luminance of the “on” period. All subjects also underwent two resting-state scans (one for perfusion and one for BOLD), during which the subject was presented with the “off” condition for 3 minutes.

Imaging Protocol

Imaging data were collected on a GE Signa Excite 3 Tesla whole body system with a body transmit coil and an eight channel receive coil. During the perfusion resting-state scan and the block design scan, data were acquired with a PICORE QUIPPS II (Wong et al. 1998) arterial spin labeling (ASL) sequence (TR=2s , TI1/TI2=600/1500ms, 10 cm tag thickness, and a 1 cm tag-slice gap) with a dual echo spiral readout (TE1/TE2=9.1/30 ms, FOV=24 cm, 64×64 matrix, and a flip angle= 90 degrees). Small bipolar crusher gradients were included to reduce signal from large vessels (b=2 s/mm2). Three oblique axial 8 mm slices were prescribed about the calcarine sulcus for this ASL run. The first echo data was used for analysis of the perfusion response and is referred to as the ASL or perfusion signal, while the second echo data was used for analysis of the BOLD response and is referred to as the 2nd echo BOLD signal. During the periodic single trial runs, BOLD-weighted images were acquired with a spiral readout (TE=25ms, TR=500 ms, FOV=24cm, 64×64 matrix, and a flip-angle of 45 degrees) and the same slice prescription as the ASL runs. The second resting state scan was acquired with the following BOLD imaging parameters (TE=25ms, TR=250 ms, FOV=24cm, 64×64 matrix, and a flip-angle of 40 degrees).

A high resolution structural scan was acquired with a magnetization prepared 3D fast spoiled gradient acquisition in the steady state (FSPGR) sequence (TI 450ms, TR 7.9ms, TE 3.1ms, 12 degree flip angle, FOV 25×25×16 cm, matrix 256×256×124).

Cardiac pulse and respiratory effort data were monitored using a pulse oximeter (InVivo) and a respiratory effort transducer (BIOPAC), respectively. The pulse oximeter was placed on the subject's left index finger, and the respiratory effort belt was placed around the subject's abdomen. Physiological data were sampled at 40 samples per second using a multi-channel data acquisition board (National Instruments).

Data Analysis

Anatomical Definition of Noise ROI

Anatomical data were segmented into gray matter, white matter, and CSF partial volume maps using the FAST algorithm available in the FSL software package (Smith et al. 2004). Tissue partial volume maps were linearly interpolated to the resolution of the functional data series using AFNI (Cox 1996). In order to form white matter ROIs, the white matter partial volume maps were thresholded at a partial volume fraction of 0.99 and then eroded by two voxels in each direction to further minimize partial voluming with gray matter. CSF voxels were determined by first thresholding the CSF partial volume maps at 0.99 and then applying a 3-dimensional nearest-neighbor criteria to minimize multiple tissue partial voluming. Since CSF regions are typically small compared to white matter regions mask, erosion was not applied.

CSF and white matter ROIs were combined to form the anatomical noise ROI. Figure 2 depicts white matter and CSF ROIs for Subject 1, as denoted by the magenta voxels, overlaid on their respective partial volume maps. We refer to the application of CompCor with the anatomical noise ROI as aCompCor.

Figure 2
Areas with a high fraction of white matter and cerebrospinal fluid (CSF), as denoted by the magenta voxels, overlaid on their respective partial volume maps from a representative slice from Subject 1. White matter-only areas (panel a) were determined ...

tSTD Based Determination of Noise ROI

In a preliminary abstract, Lund and Hanson (2001) showed that voxel time courses with a relatively high temporal standard deviation were dominated by physiological noise. They observed that these voxels occurred in edge-regions, ventricles, and in areas close to large vessels. In their approach, they manually selected five pixels with high temporal standard deviation (tSTD) that appeared to represent physiological fluctuation. The time-series from these voxels were then included as nuisance covariates in a GLM, resulting in a marked improvement in detection power (Lund and Hanson 2001). Here we extend the prior work by first using the temporal standard deviation to select voxels in an unsupervised fashion and then using principal components analysis to reduce the dimensionality of the data. For each voxel time series, the temporal standard deviation is defined as the standard deviation of the time series after the removal of low frequency nuisance terms (e.g. linear and quadratic drift). We denote the noise ROI determined in this fashion as the tSTD noise ROI and refer to the application of CompCor with this ROI as tCompCor.

To confirm that the tSTD metric identifies voxels with high levels of physiological noise, we first examined the relation between measures of physiological noise and tSTD. For each voxel, we defined the fractional variance of physiological noise as the ratio of the variance of the voxel time series accounted for by physiological noise regressors as determined with RETROICOR (see section on GLM analysis below) to the variance of the original time series after removal of constant and linear trends (e.g. the square of tSTD). The fractional variance of physiological noise was then compared to the temporal standard deviation on a per-voxel basis. Figure 3 shows spatial maps of the fractional variance and tSTD (panels a and b respectively) as well as a plot of tSTD versus fractional variance for a representative slice from subject 1. As demonstrated by the plot in panel c and the similarity of the spatial maps, the voxels with the highest tSTD also tend to have the highest fractional variance of physiological noise, confirming the observations of (Lund and Hanson 2001).

Figure 3
Panel (a) shows a spatial map of the fractional variance of physiological noise for the resting BOLD scan from a representative slice in subject 1. Panel (b) shows a spatial map of the temporal standard deviation (tSTD). Areas of high fractional variance ...

To construct the tSTD noise ROI, we sorted the voxels by their temporal standard deviation and retained a pre-specified upper fraction of the sorted voxels within each slice. For example, specification of an upper fraction of 1% retains the voxels with the top 1% of temporal standard deviation values. Specification of the fraction involves a trade-off between including too few voxels to accurately represent the physiological noise space versus including too many voxels, which would tend to include voxels whose temporal standard deviation is not dominated by physiological noise. To empirically determine a reasonable value to use for the upper fraction, we computed the mean fractional variance of physiological noise across voxels within the tSTD ROI as a function of the upper fraction. This analysis was performed for each subject for both resting ASL and BOLD runs. Figure 4 presents the mean across subjects of the fractional variance of physiological noise versus the fraction of included voxels. As a greater number of voxels are included in the analysis, the average fractional variance of physiological noise decreases for both ASL and BOLD resting runs. Both the ASL and BOLD curves show a steady decline in the fractional variance above a threshold of approximately 2%. Based on these results, we chose a 2% threshold (~20−30 voxels per slice) as a reasonable empirical threshold that effectively identified voxels with the highest fractional variance of physiological noise.

Figure 4
Mean fractional variance of physiological noise across subjects as a function of the fraction of voxels (sorted by tSTD) that is included in the noise ROI for BOLD (red) and ASL (blue) resting-state data. The black dotted line denotes the 2% threshold ...

Exclusion of stimulus-related components

Although the construction of the noise ROIs is designed to include voxels that are unlikely to include stimulus-correlated activity, there is always a finite probability that some voxels may contain stimulus-related components. Inclusion of these voxels in the noise ROI will tend to reduce the statistical performance of the CompCor algorithm because stimulus-related components will be treated as physiological noise components. To minimize the probability of including voxels with stimulus-related components, we correlated the time course of each voxel within the noise ROI with the stimulus-related reference function and excluded any voxels with a p-value less than 0.2. To assess the effect of this threshold on performance, we performed Monte Carlo simulations to generate receiver operating characteristic (ROC) curves at various threshold levels and used the area under the ROC curve to quantify performance. An area of 1 represents ideal performance.

Simulated voxel responses (N=5000) were generated using the BOLD GLM presented in equation 1. Physiological regressors, P, were generated at frequencies consistent with cardiac (0.9 Hz) and respiration (0.3 Hz) with uniformly distributed phases to represent phase lags between voxels. Physiological noise weights, c, were chosen from a normal distribution N(0, 0.3). A constant term was included as a nuisance term with the regressor weight d chosen from a uniform distribution U(0,1), and the additive white noise term, n, was generated from a normal distribution N(0,1). The block stimulus paradigm was used to construct the stimulus design matrix X, and the stimulus weight h was set to a value of either 0 or 0.3 for the null and alternative hypotheses, respectively. Correlation of the simulated voxel responses with the block stimulus response was used to calculate p-values and z-scores for each condition. In the null hypothesis condition, the distribution of z-scores had a mean and standard deviation of −0.31 and 0.89, respectively, whereas in the alternative hypothesis condition the mean and standard deviation were 3.18 and 1.27 , respectively.

To construct the ROC curve, we varied the threshold on p-values from 0.01 to 1.0. For each threshold value, voxel time-series with a p-value above the threshold were placed column-wise into a matrix M. A principal component analysis of the matrix M (see details in following section) was then used to form an estimate of the physiological noise matrix Pest. A GLM analysis of the simulated voxel time-series was performed using the idealized stimulus response X and the computed Pest. The calculated p-values obtained under the null and alternative hypothesis conditions were then used to generate a ROC curve for each threshold value.

As the exclusion threshold was reduced from p = 1.0 down to p = 0.1, there was little degradation in detector performance (area under the ROC curve >0.99), since voxels with any hint of activation were effectively excluded. As the threshold was further reduced from 0.1 to 0.05, the ROC area decreased from 0.99 down to a value of 0.96, which is less than the area (0.98) obtained in the absence of CompCor (e.g. Pest not included in the GLM). The degradation in performance with lower p-values reflects the inclusion of weak stimulus-correlated components. Based on these findings, we chose a conservative threshold of p = 0.2 to exclude voxels with stimulus-correlated fluctuations.

Determination of Principal Components

Voxel time-series from the noise ROI (either anatomical or tSTD) were placed in a matrix M of size N × m, with time along the row dimension and voxels along the column dimension. The constant and linear trends of the columns in the matrix M were removed prior to column-wise variance normalization. The covariance matrix C = MMT was constructed and decomposed into its principal components using a singular value decomposition.

The number of significant principal components to retain was determined using a modified version of the “broken stick” method described in (Jackson 1993). This method identifies meaningful principal components (e.g. components unlikely to be due to random noise) by comparing their associated principal values to principal values derived from normally distributed data of equal rank. Ordered principal values calculated from normally distributed data tend to show a sharp initial decrease followed by a more gradual decrease, thus resembling a “broken-stick”. In this method, a Monte-Carlo simulation (N = 1000) was first used to generate a statistical representation of expected principal values derived from normally distributed data of rank equal to the matrix M. This statistical distribution was then compared to the computed principal values from the data. Principal components that were significantly larger than the generated distribution, as assessed using a two-tailed t-test (p<0.05), were retained. Based on this method, we found an average of 6.3 ± 0.52 and 4.5 ± 0.38 significant principal components for BOLD and ASL runs, respectively, when using the anatomical noise ROI. For the tSTD noise ROI there were an average of 5.9 ± 0.74 and 4.2 ± 0.59 principal components for the BOLD and ASL runs, respectively.

General Linear Model Analysis

All images were coregistered using AFNI software (Cox 1996). Statistical inference was performed using a general linear model (GLM) analysis with an AR(1) model for the additive noise component (Burock and Dale 2000; Woolrich et al. 2001). The analysis was performed using the appropriate GLM for the BOLD and ASL time-series under the following conditions: 1) no inclusion of physiological noise regressors (e.g. no correction), 2) inclusion of RETROICOR regressors, 3) inclusion of CompCor regressors derived from the anatomical noise ROI (aCompCor), and 4) inclusion of CompCor regressors derived from the tSTD noise ROI (tCompCor). The stimulus-related component Xh was formed by convolving the appropriate stimulus pattern with a gamma density function of the form h(t) = (τn!)− 1((t − Δt)/τ) exp(− (t − Δt)/τ) for t ≥ Δt and 0 otherwise, with τ = 1.2, n = 3 and Δt = 1 (Boynton et al. 1996). As nuisance regressors, we used a constant term, a linear term, and a discrete cosine set implemented in SPM2 with a minimum period of 120 seconds (Worsley and Friston 1995; Lund et al. 2006). For statistical inference with RETROICOR, physiological noise regressors obtained from the cardiac and respiratory measurements were used to form the physiological noise matrix P (Glover et al. 2000b; Restom et al. 2006). In addition, estimates of the cardiac and respiratory components were formed from partitioning the estimate of Pc (i.e., PCcC and PRcR where the C and R subscripts denote cardiac and respiratory, respectively). For the application of aCompCor and tCompCor, GLM analysis was performed using noise matrix Pest, constructed from the principal component analysis of voxel time-series from the anatomical noise ROI and the tSTD noise ROI, respectively. The estimate of the term Pestc was used as the CompCor estimate of the physiological noise component.

For each method, we formed a “corrected time series” by subtracting the estimates of the physiological noise components and the nuisance parameters from the measured time series. We then used these corrected time series to compute, on a per-subject basis, the mean temporal standard deviation of the resting-state data across all voxels within gray matter (partial volume > 0.9). Paired t-tests (two-tailed) were used to assess differences between methods across the sample of subjects.

The number of activated voxels detected with each method was computed for each subject by thresholding the F-statistics at 5 and 30 for the ASL and BOLD functional runs, respectively. These values were chosen to yield approximately the same number of activated voxels for the ASL and BOLD functional runs and are consistent with thresholds previously used to investigate physiological noise reduction for ASL (Restom et al. 2006). Paired t-tests (two-tailed) were used to assess differences in the number of activated voxels.

Spectral Analysis

We used spectral analysis to compare the physiological components estimated by RETROICOR with those estimated by CompCor. For each subject, we computed the power spectra of the physiological noise components estimated with each method on a per-voxel basis, normalized each power spectrum by its peak value, and then averaged the normalized spectra across all voxels in gray matter (partial volume > 0.9). The RETROICOR power spectra were used to identify the peak cardiac and respiratory frequencies. We defined cardiac and respiratory frequency bands as the 0.1 Hz wide band of frequencies centered around the respective peak frequency. To quantify the similarity between methods of the spectra in these bands, we computed the band-averaged coherence Cohxy (fband) defined as

equation M3

where fband denotes the frequency band, x and y are the two time series of interest, Sxy(f) is the cross-spectrum of x and y, and Sxx(f) and Syy(f) are the power spectra of x and y, respectively (Sun et al. 2004). The cross-spectral and power-spectral densities are estimated using Welch's modified periodogram averaging method as previously outlined by (Sun et al. 2004). Coherence is bounded by 0 and 1, with 1 representing perfect coherence, and is analogous to correlation analysis of time-series. The spectral coherence between physiological noise elements was computed on a per-voxel basis and then averaged across gray matter voxels (partial volume > 0.9).

Additional Analysis of BOLD Data with TR of 2 Seconds

Although the physiological noise in the BOLD time series is critically sampled given the short TR of 0.5s, the cardiac and respiratory noise fluctuations in the ASL perfusion time series will typically be aliased due to the lower sampling rate (TR = 2s). In the results section we show that both aCompCor and tCompCor are capable of identifying these aliased components. To further demonstrate the performance of aCompCor and tCompCor for undersampled data, we downsampled the BOLD time series by a factor of 4 to generate a BOLD time series with a TR of 2 seconds and analyzed the performance of the CompCor methods using the GLM analysis approach described above. We also applied CompCor to the BOLD weighted second echo data from the ASL time series (TR of 2 s), using the GLM in Equation 2 with hBOLD as the parameter of interest.

Receiver Operating Characteristic Curve

To gain additional insight into the relative performance of the CompCor methods, we used simulations to construct receiver operating characteristic (ROC) curves. For these simulations, we first defined a region of interest (ROI) in the visual cortex based on GLM analysis (with tCompCor regressors) of the 2nd echo BOLD data from the block design functional runs. For each subject, a liberal threshold (r>0.2) was used to delimit an ROI with approximately 400 to 500 voxels. These voxels were then used for the subsequent simulations, which made use of the BOLD resting-state data. For each method, the false positive rates at varying thresholds were calculated by applying the GLM analysis to the resting-state data and computing the fraction of false positives within the ROI. To determine the true positive rates, we added simulated activations (block design: 20s on, 40 seconds off; 1% amplitude) to the time course of each voxel within the ROI. At each threshold level, we then computed the fraction of true positives within the ROI.

Performance of tCompCor in the Presence of Motion

As the tCompCor approach is based on the identificaton of voxel time series with high temporal standard deviations, it will also identify voxels with large signal components due to subject motion. Motion-related signal changes can reflect both simple linear effects (e.g. translation and rotation of brain regions) and more complex nonlinear effects, such as changes in the distortion and blurring of the image due to magnetic field inhomogeneities. To the extent that motion-related signal changes increase the temporal standard deviation of the time series, these changes can be reflected in the principal components identified by tCompCor. In the presence of stimulus-correlated motion, the exclusion process described previously will bias against the inclusion of time series with stimulus-correlated changes. As a demonstration of the application of tCompCor in the presence of subject motion, we analyzed a BOLD (TR = 0.5s) time series in which there was significant subject motion. The data were acquired during a periodic single trial run with the protocol and imaging parameters described above. To emphasize the effect of motion, we applied tCompCor to data that were not motion corrrected. Estimates of rotation and displacement during the experiment were obtained with the 3dvolreg program in AFNI.

Results

Figure 5 shows normalized power spectra of physiological components estimated by RETROICOR, aCompCor, and tCompCor for the resting-state BOLD run of subject 1. The top row (panel a) depicts the respiratory and cardiac components identified with the use of RETROICOR. Cardiac and respiratory peaks are located prominently at ~1.2 Hz and ~0.2 Hz, respectively. Panels b) and c) show the average spectra of elements estimated by aCompCor and tCompCor, respectively. Both variants of CompCor estimate cardiac and respiratory elements that are similar to those removed by RETROICOR. In addition, very low frequency 1/f components are visible in both the aCompCor and tCompCor spectra. These low frequency components can appear in the CompCor spectra because we remove only constant and linear trend terms prior to the principal components analysis. Note that, as stated in the Methods section, the use of the discrete cosine transform (DCT) low frequency nuisance terms is reserved for the statistical assessments performed with the GLM. As a result, the space spanned by the principal components may partially overlap with the space spanned by the discrete cosine transform nuisance terms. Potential issues with this overlap are addressed in the Discussion section.

Figure 5
Average normalized power spectra of components estimated with the application of various correction schemes to the resting BOLD run from subject 1. As shown in panel (a), cardiac and respiratory elements estimated by RETROICOR are located at 1.2 and 0.2 ...

Figure 6 shows the average spectra from Subject 1 of the noise elements estimated from the resting-state ASL data. In the top row, the RETROICOR estimated cardiac (red line) and respiratory (green line) elements are aliased due to the long TR. The spectrum of elements estimated by either aCompCor (panel b) and tCompCor (panel c) show components consistent with the cardiac and respiratory elements identified by RETROICOR.

Figure 6
Average normalized power spectra of components estimated by the various correction schemes for the resting ASL run from subject 1. As shown in panel (a), cardiac (red) and respiratory (green) elements identified by RETROICOR are aliased due to the long ...

To quantify the similarity between the spectra, we used the spectral coherence analysis described above in Methods. For aCompCor applied to the resting BOLD runs, the average spectral coherence and standard error values across subjects (N=10) were 0.67 ± 0.08 and 0.80 ± 0.05 for respiratory and cardiac elements, respectively. Corresponding values for tCompCor were 0.62 ± 0.10 and 0.86 ± 0.07. For resting ASL runs the mean spectral coherences when using aCompCor were 0.79 ± 0.03 and 0.65 ± 0.06 for respiratory and cardiac elements, respectively, with corresponding values of 0.82 ± 0.03 and 0.69 ± 0.06 for tCompCor.

The effects of CompCor and RETROICOR on the average temporal standard deviation (tSTD) of the resting BOLD and ASL data are shown in Figure 7, with associated standard error bars. The panels show the normalized mean temporal standard deviation (tSTD) in gray matter (partial voluming >0.9) for a) resting TR = 0.5 s BOLD data, b) downsampled (TR = 2s) resting BOLD data, c) resting 1st echo ASL (perfusion) data, and d) resting 2nd echo BOLD data, with the application of RETROICOR and the two variants of CompCor. For each subject, the tSTD values obtained with physiological noise reduction methods are normalized by the values obtained for the uncorrected time series. The standard error bars are provided to give a sense of the inter-subject variability, but should not be used to interpret statistical significance, since paired t-tests are used to compare the normalized tSTD values between methods.

Figure 7
Percent temporal standard deviation across gray matter voxels (partial volume >0.9) for uncorrected data (denoted as None) and data after application of RETROICOR (denoted as Phys) and CompCor for (a) resting-state BOLD data, (b) downsampled resting-state ...

As compared to no correction, RETROICOR significantly (p<0.001) reduced the normalized tSTD in the resting (TR = 0.5s) BOLD data by 8% whereas application of aCompCor and tCompCor resulted in even greater reductions of 20% (p<0.001) and 29% (p<0.001), respectively. Both variants of CompCor significantly (p<0.02) reduced the tSTD as compared to RETROICOR, and tCompCor significantly (p<0.001) reduced the tSTD compared to aCompCor. For the downsampled (TR = 2s) BOLD data, all three methods significantly (p< 0.001) reduced the normalized tSTD with percent decreases of 12% , 17% and 22% for RETROICOR, aCompCor, and tCompCor, respectively. In addition, the normalized tSTD values achieved with tCompCor were significantly (p<0.03) lower than the values obtained with either RETROICOR or aCompCor.

For the resting 1st echo ASL data, all three methods significantly (p< 0.001) reduced the normalized tSTD with percent decreases of 13%, 18%, and 23%. The normalized tSTD values obtained with aCompCor and tCompCor were significantly (p<0.04) lower than the values obtained with RETROICOR. In addition, the tSTD obtained with tCompCor were significantly (p< 0.03) lower than the values obtained with aCompCor. Similarly, for the the resting 2nd echo BOLD data, the three methods led to significant (p< 0.001) reductions in tSTD with percent decreases of 13%, 19%, and 21%. Both aCompCor and tCompCor lead to significantly greater reductions (p< 0.01) than RETROICOR.

Bar graphs comparing the effect of the various correction schemes on the normalized number of activated voxels for the BOLD periodic and functional ASL block runs are shown in Figure 8. For each method, the number of activated voxels is normalized on a per subject basis by the number of voxels detected for the data prior to the application of noise reduction. As with the plots in Figure 7, since statistical signifiance is assessed with paired tests, the standard errors bars are provided only to give a sense of inter-subject variability. As compared to no correction, there were significant (p< 0.02) increases of 10%, 76% and 82% in the normalized number of activated voxels for the periodic (TR = 0.5s) BOLD runs (panel a) with the application of RETROICOR, aCompCor, and tCompCor, respectively. Application of either aCompCor or tCompCor increased the number of voxels significantly (p< 0.04) as compared to RETROICOR. For the downsampled (TR = 2s) BOLD data, both aCompCor and tCompCor significantly increased the normalized number of activated voxels with respect to both the uncorrected data (p< 0.02) and RETROICOR (p< 0.02). In addition, the number of activated voxels with tCompCor was significantly greater (p< 0.04) than the number obtained with aCompCor.

Figure 8
Percent of significantly activated voxels across subjects (N=10) for uncorrected data (denoted as None) and data with application of RETROICOR (denoted as Phys) and CompCor for (a) periodic design BOLD data, (b) downsampled periodic design BOLD data, ...

For the first echo ASL block functional data, the normalized number of activated voxels increased significantly (p<0.03) with gains of 65%, 51%, and 49% with the application of RETROICOR, aCompCor, and tCompCor, respectively. There was not a significant difference (p> 0.25) between the normalized number of activated voxels detected with the three methods. For the second echo BOLD functional data, both aCompCor and tCompCor significantly increased the normalized number of activated voxels with respect to both the uncorrected data (p< 0.005) and RETROICOR (p< 0.005).

Figure 9 shows a representative ROC curve comparing the relative performances achieved without correction (blue) and the application of RETROICOR (green), aCompCor (red), and tCompCor (cyan). Consistent with the results shown in Figure 8, both aCompCor and tCompCor yield better detection performance than RETROICOR, and the performance of tCompCor is better than that of aCompCor.

Figure 9
Representative receiver operating characteristic curve showing the true positive rate versus false positive rate for uncorrected data (blue), RETROICOR (green), aCompCor (red), and tCompCor (cyan).

Figure 10 shows an example of the performance of tCompCor in the presence of signal changes due to subject motion. The top row shows the uncorrected BOLD time series (red), the time series after application of tCompCor (blue), and the periodic single trial reference function (black dash line). For this run, the most significant motion components were found to be left-right roll and displacement. Estimates of roll in units of degrees (blue) and displacement in units of millimeters (green) are shown in the middle row. Note that the large motion peaking at 79 seconds gives rise to a large signal change in the uncorrected time series. The bottom row shows the top three principal components as estimated by tCompCor. These components reflect the motion at 79 seconds as well as smaller movements in other portions of the experimental run. Projecting these components out of the uncorrected time series leads to the corrected time series, with a noticeable reduction in the motion-related signal changes, especially at the 79 second timepoint.

Figure 10
Example of the application of tCompCor in the presence of motion-related signal changes. The top row shows the original time series (red), the time series after the application of tCompCor (blue), and the periodic design reference function (black dashed). ...

Discussion

In this paper, we have examined whether signal components derived from regions of interest that are unlikely to be modulated by neural activity can be used to estimate noise components (due to physiological fluctuations, subject motion, etc.) within activated regions. We considered two methods for the determination of the noise ROIs: 1) anatomical identification of significant areas of CSF and white matter and 2) definition of noise regions based upon their temporal standard deviation. We demonstrated that the application of CompCor using either ROI significantly reduces the temporal standard deviation of both resting state BOLD and ASL data. Additionally, we have shown that CompCor using either ROI leads to a marked improvement in sensitivity to detect the response to a visual stimulus, as quantified by the number of significantly activated voxels for both BOLD and ASL experiments.

For the resting-state data, we found a significantly greater reduction of the normalized temporal standard deviation with the application of CompCor as compared to RETROICOR (with the exception of aCompCor for the downsampled BOLD data). This most likely due to the removal of noise terms that are identified by the principal components analysis used in CompCor but not modeled by the RETROICOR cardiac and respiratory regressors.

For the functional BOLD runs, we found that the application of either aCompCor or tCompCor resulted in significantly more activated voxels as compared to RETROICOR. In contrast, for the ASL data, there were not significant differences between the percent of activated voxels detected by either form of CompCor as compared to RETROICOR. Due to the tag and control modulation used in the ASL acquisition process, low frequency noise components are effectively attenuated in the ASL analysis process (Liu and Wong 2005). The fact that the improved detection performance of the CompCor methods was limited to the BOLD data suggests that the gain reflects a reduction in low frequency noise components that are identified by CompCor but not by RETROICOR. As noted in the Results section, the presence of these low frequency components may cause the space spanned by the principal components to overlap with the signal space spanned by the nuisance regressors used in the general linear model. In the unlikely case that the overlap causes the model to be ill-posed, it may be necessary to orthogonalize the principal components with respect to the nuisance regressors. However, we have not found this to be necessary in practice.

The performance gains achieved with both RETROICOR and CompCor (as compared to the uncorrected data) were observed both for the critically sampled (TR = 0.5s) BOLD data and for the undersampled or downsampled (TR = 2s) ASL and BOLD data. In the undersampled data, both respiratory and cardiac fluctuations are typically aliased. With RETROICOR, these aliased components can be identified from the external measures of cardiac and repiratory activity, as has been previously discussed (Glover et al. 2000b; Restom et al. 2006). In CompCor, the principal components analysis takes advantage of the spatial coherence of the aliased components to identify the most significant aliased components. Because this identification procedure does not depend on the frequencies of the physiological components, it works well even with undersampled data. As shown by Figure 6 and the spatial coherence measures reported in the Results section, both RETROICOR and CompCor identify similar spectral components in both the critically sampled and undersampled data.

The ability of CompCor to identify spatially coherent noise components that are unlikely to be of neural origin also plays a role in its capacity to improve performance in the presence of motion-related signal changes. Since these effects are primarily due to bulk motion, they will tend to give rise to spatially coherent signal changes. Although the exact signal change in any given voxel will depend on a variety of factors, such as the local anatomy and magnetic field distributions, CompCor will identify those components that are common across voxels. The performance of CompCor in any given voxel will depend on the extent to which these components can represent the signal changes in that voxel. An example of the efficacy of tCompCor in removing a motion-related signal component is shown in Figure 10. In cases with especially severe motion artifacts, the principal components identified by CompCor may primarily describe the effects of motion, thus possibly decreasing the ability of CompCor to reduce the effect of cardiac and repiratory fluctuations. Further work to better elucidate the limitations of CompCor in the presence of motion would be useful.

The accurate specification of the noise ROI is key to both aCompCor and tCompCor. As discussed in the methods section, inclusion of voxels with stimulus-correlated activity can reduce the performance of the methods. For aCompCor this may occur due to factors such as inaccuracies in tissue segmentation or misalignment of the anatomical and functional data (e.g. due to subject movement between the scans). For tCompCor this performance reduction can occur when voxels with high tSTD also have significant stimulus-correlated components (e.g. due to stimulus-correlated motion). In the present study, we used a liberal threshold (p< 0.2) to exclude voxels with possible stimulus-correlated components. Although our experimental data support the efficacy of this approach for periodic and block designs, the performance of this approach may be reduced for complex event-related fMRI experiments in which the expected stimulus-related response may not be well defined. In addition, in the rare cases where the cardiac noise is aliased to the stimulus frequency (Lund and Hanson 2001), the exclusion criteria will reduce the ability of CompCor to identify these noise sources. Finally, this exclusion approach is not easily extended to resting state experiments aimed at studying functional connectivity. Future work addressing the applicability of CompCor to these types of studies would be useful.

In this work we used an empirical threshold to determine the fraction of voxels to retain in the tSTD noise ROI. While this threshold provided better performance for the data sets analyzed in this study, it is possible that this threshold may not be optimum for other data sets. An examination of additional datasets would be helpful for better characterizing the efficacy of tCompCor. In addition, an investigation into other approaches for determining the threshold may provide gains in performance.

Conclusion

We have shown that application of CompCor to ASL and BOLD fMRI time-series can significantly reduce noise due to physiological fluctuations and other sources, such as subject motion. CompCor does not require external monitoring and can be applied in an automated fashion to reduce the confounding effect of physiological fluctuations on fMRI time-series.

Acknowledgments

The authors thank Peter Costandi and Joanna Perthen for their valued assistance with the preparation of this paper. This work was supported in part by a Biomedical Engineering Research grant from the Whitaker Foundation and by NIH Grant R01NS051661.

Footnotes

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.

References

  • Beckmann CF, Smith SM. Probabilistic independent component analysis for functional magnetic resonance imaging. IEEE Trans Med Imaging. 2004;23(2):137–52. [PubMed]
  • Birn RM, Diamond JB, Smith MA, Bandettini PA. Separating respiratory-variation-related fluctuations from neuronal-activity-related fluctuations in fMRI. Neuroimage. 2006;31(4):1536–48. [PubMed]
  • Biswal B, DeYoe AE, Hyde JS. Reduction of physiological fluctuations in fMRI using digital filters. Magn Reson Med. 1996;35(1):107–13. [PubMed]
  • Boynton GM, Engel SA, Glover GH, Heeger DJ. Linear systems analysis of functional magnetic resonance imaging in human V1. J. Neuroscience. 1996;16:4207–4221. [PubMed]
  • Burock MA, Dale AM. Estimation and detection of event-related fMRI signals with temporally correlated noise: a statistically efficient and unbiased approach. Human Brain Mapping. 2000;11:249–260. [PubMed]
  • Cox RW. AFNI-software for analysis and visualization of functional magnetic resonance neuroimages. Comput. Biomed. Res. 1996;29:162–173. [PubMed]
  • Dagli MS, Ingeholm JE, Haxby JV. Localization of cardiac-induced signal change in fMRI. Neuroimage. 1999;9(4):407–15. [PubMed]
  • Glover GH, Li T-Q, Ress D. Image-based method for retrospective correction of physiological motion effects in fMRI: RETROICOR. Magn Res Med. 2000a;44:162–167. [PubMed]
  • Glover GH, Li TQ, Ress D. Image-based method for retrospective correction of physiological motion effects in fMRI: RETROICOR. Magn Reson Med. 2000b;44(1):162–7. [PubMed]
  • Hoge RD, Atkinson J, Gill B, Crelier GR, Marrett S, Pike GB. Investigation of BOLD signal dependence on cerebral blood flow and oxygen consumption: the deoxyhemoglobin dilution model. Magn Reson Med. 1999;42(5):849–63. [PubMed]
  • Hu X, Le TH, Parrish T, Erhard P. Retrospective estimation and correction of physiological fluctuation in functional MRI. Magn Reson Med. 1995;34(2):201–12. [PubMed]
  • Jackson DA. Stopping Rules in Principal Components Analysis: A Comparison of Heuristical and Statistical Approaches. Ecology. 1993;74(8):2204–2214.
  • Josephs O, Howseman A, Friston KJ, Turner R. Physiological Noise Modelling for multi-slice EPI fMRI using SPM. Proc. Intl. Soc. Mag. Reson. Med. 2001:1682.
  • Kruger G, Glover GH. Physiological noise in oxygenation-sensitive magnetic resonance imaging. Magn Reson Med. 2001;46(4):631–7. [PubMed]
  • Liu CS, Miki A, Hulvershorn J, Bloy L, Gualtieri EE, Liu GT, Leigh JS, Haselgrove JC, Elliott MA. Spatial and temporal characteristics of physiological noise in fMRI at 3T. Acad Radiol. 2006;13(3):313–23. [PubMed]
  • Liu TT, Wong EC. A signal processing model for arterial spin labeling functional MRI. Neuroimage. 2005;24(1):207–15. [PubMed]
  • Lund TE, Hanson LG. Physiological Noise Reduction in fMRI Using Vessel Time-Series as Covariates in a General Linear Model. Proc. Intl. Soc. Mag. Reson. Med. 2001;9:22.
  • Lund TE, Madsen KH, Sidaros K, Luo WL, Nichols TE. Non-white noise in fMRI: does modelling have an impact? Neuroimage. 2006;29(1):54–66. [PubMed]
  • McKeown MJ, Hansen LK, Sejnowsk TJ. Independent component analysis of functional MRI: what is signal and what is noise? Curr Opin Neurobiol. 2003;13(5):620–9. [PMC free article] [PubMed]
  • Mumford JA, Hernandez-Garcia L, Lee GR, Nichols TE. Estimation efficiency and statistical power in arterial spin labeling fMRI. Neuroimage. 2006;33(1):103–14. [PMC free article] [PubMed]
  • Petersen N, Jensen J, Burchardt J, Stodkilde-Jorgensen H. State space models for physiological noise in fMRI time series. Neuroimage. 1998;7(4):727–741.
  • Pfeuffer J, Van de Moortele PF, Ugurbil K, Hu X, Glover GH. Correction of physiologically induced global off-resonance effects in dynamic echo-planar and spiral functional imaging. Magn Reson Med. 2002;47(2):344–53. [PubMed]
  • Restom K, Behzadi Y, Liu TT. Physiological noise reduction for arterial spin labeling functional MRI. Neuroimage. 2006;31(3):1104–15. [PubMed]
  • Smith SM, et al. Advances in functional and structural MR image analysis and implementation as FSL. Neuroimage. 2004;23(Suppl 1):S208–19. [PubMed]
  • Sun FT, Miller LM, D'Esposito M. Measuring interregional functional connectivity using coherence and partial coherence analyses of fMRI data. Neuroimage. 2004;21(2):647–58. [PubMed]
  • Thomas CG, Harshman RA, Menon RS. Noise reduction in BOLD-based fMRI using component analysis. Neuroimage. 2002;17(3):1521–37. [PubMed]
  • Wong EC, Buxton RB, Frank LR. Quantitative imaging of perfusion using a single subtraction (QUIPSS and QUIPSS II). Magn. Reson. Med. 1998;39(5):702–708. [PubMed]
  • Woolrich MW, Ripley BD, Brady M, Smith SM. Temporal autocorrelation in univariate linear modeling of FMRI data. Neuroimage. 2001;14(6):1370–86. [PubMed]
  • Worsley KJ, Friston KJ. Analysis of fMRI time-series revisted --- again. Neuroimage. 1995;2:173–181. [PubMed]