Differentiating BOLD and Non-BOLD Signals in fMRI Time Series Using Multi-Echo EPI (original) (raw)

Neuroimage. Author manuscript; available in PMC 2013 Apr 15.

Published in final edited form as:

PMCID: PMC3350785

NIHMSID: NIHMS363220

Prantik Kundu

§Section on Functional Imaging Methods, Laboratory of Brain and Cognition, National Institute of Mental Health, National Institutes of Health, Bethesda, MD, USA 20892

£Department of Psychiatry, University of Cambridge, Addenbrooke’s Hospital, Hills Road, Cambridge, UK CB2 2QQ

Souheil J. Inati

¶Functional MRI Facility, National Institute of Mental Health, National Institutes of Health, Bethesda, MD, USA 20892

Jennifer W. Evans

§Section on Functional Imaging Methods, Laboratory of Brain and Cognition, National Institute of Mental Health, National Institutes of Health, Bethesda, MD, USA 20892

⋄Center for Neuroscience and Regenerative Medicine, Henry M. Jackson Foundation, Rockville, Maryland, USA 20852

Wen-Ming Luh

¶Functional MRI Facility, National Institute of Mental Health, National Institutes of Health, Bethesda, MD, USA 20892

Peter A. Bandettini

§Section on Functional Imaging Methods, Laboratory of Brain and Cognition, National Institute of Mental Health, National Institutes of Health, Bethesda, MD, USA 20892

¶Functional MRI Facility, National Institute of Mental Health, National Institutes of Health, Bethesda, MD, USA 20892

§Section on Functional Imaging Methods, Laboratory of Brain and Cognition, National Institute of Mental Health, National Institutes of Health, Bethesda, MD, USA 20892

¶Functional MRI Facility, National Institute of Mental Health, National Institutes of Health, Bethesda, MD, USA 20892

£Department of Psychiatry, University of Cambridge, Addenbrooke’s Hospital, Hills Road, Cambridge, UK CB2 2QQ

⋄Center for Neuroscience and Regenerative Medicine, Henry M. Jackson Foundation, Rockville, Maryland, USA 20852

Corresponding Author: Prantik Kundu, Tel: 301-402-1342, Fax: 301-402-1370, vog.hin.liam@pudnuk

Abstract

A central challenge in the fMRI based study of functional connectivity is distinguishing neuronally related signal fluctuations from the effects of motion, physiology, and other nuisance sources. Conventional techniques for removing nuisance effects include modeling of noise time courses based on external measurements followed by temporal filtering. These techniques have limited effectiveness. Previous studies have shown using multi-echo fMRI that neuronally related fluctuations are Blood Oxygen Level Dependent (BOLD) signals that can be characterized in terms of changes in R2* and initial signal intensity (S0) based on the analysis of echo-time (TE) dependence. We hypothesized that if TE-dependence could be used to differentiate BOLD and non-BOLD signals, non-BOLD signal could be removed to denoise data without conventional noise modeling. To test this hypothesis, whole brain multi-echo data were acquired at 3 TEs and decomposed with Independent Components Analysis (ICA) after spatially concatenating data across space and TE. Components were analyzed for the degree to which their signal changes fit models for R2* and S0 change, and summary scores were developed to characterize each component as BOLD-like or not BOLD-like. These scores clearly differentiated BOLD-like “functional network” components from non BOLD-like components related to motion, pulsatility, and other nuisance effects. Using non BOLD-like component time courses as noise regressors dramatically improved seed-based correlation mapping by reducing the effects of high and low frequency non-BOLD fluctuations. A comparison with seed-based correlation mapping using conventional noise regressors demonstrated the superiority of the proposed technique for both individual and group level seed-based connectivity analysis, especially in mapping subcortical-cortical connectivity. The differentiation of BOLD and non-BOLD components based on TE-dependence was highly robust, which allowed for the identification of BOLD-like components and the removal of non BOLD-like components to be implemented as a fully automated procedure.

Keywords: fMRI, Resting state, Multi-echo, BOLD, TE dependence, ICA, Denoising, Subcortical, Connectivity

Introduction

The exploration and use of “resting state” functional magnetic resonance imaging (fMRI) data has shown explosive growth in recent years. Methods for analyzing resting state fMRI functional networks include the seed-voxel approach, which involves calculation of the correlation between a signal time course from one brain region and the time courses from the rest of the brain (Biswal et al., 1995), and the decomposition approach, involving the use of techniques such as Independent Components Analysis (ICA; Beckmann and Smith, 2004b). The consistency of resting networks in healthy adults is well established (Damoiseaux, 2006a; review by Van Den Heuvel and Pol, 2010), and the variations of networks due in several neuropsychiatric conditions have been studied (review by Broyd et al., 2009). Functional connectivity analysis of data from standard (i.e. single-echo) fMRI pulse sequences is limited, however, by the fundamental problem that in such experiments, Blood Oxygen Level Dependent (BOLD) signal arising from spontaneous neuronal activity is not differentiable from fluctuations arising from cardiac and respiratory physiology, motion, and many other sources. Several techniques have been developed to remove artifactual signal, including the use of temporal noise models and band pass filtering (Jonsson et al., 1999), but these approaches are limited in their effectiveness (Birn et al., 2006). For this reason, current de-noising techniques can underestimate the effect of non-neuronal fluctuations, remove neuronally related fluctuations, and require assumptions that may not be consistent across subjects and scan sites. Here, building on Peltier and Noll (Peltier and Noll, 2000) and Kruger and Glover (Krüger and Glover, 2001), we introduce a new method that employs multi-echo acquisition and a TE-dependence test to remove artifactual fluctuations more effectively than these previous approaches by cleanly separating BOLD and non-BOLD signal components of resting state data.

A change in BOLD contrast can be described as a change in the transverse relaxation rate or R2* due to changes in blood oxygenation (Ogawa, 1990a; Menon et al., 1993; Bandettini et al., 1994). The use of this endogenous contrast has become the primary MRI-based method by which brain activation is assessed (Ogawa, 1990b; Kwong et al., 1992; Bandettini et al., 1992). Typical acquisition of time series data involves a single TE that is equal to the resting transverse relaxation time, T2*, equal to 1/R2* (Menon et al., 1993). Over the years, multi-echo acquisition has been used to enhance our understanding of fMRI time series. Early studies acquired multi-echo fMRI during motor and visual stimulation to evaluate the effect of TE on activation mapping (Barth et al., 1999). They showed that TE influences sensitivity to various vessel effects, which expanded on findings from previous studies (Bandettini et al., 1994). Other multi-echo studies characterized baseline and activation-induced changes in R2*, which was not possible with time series acquisition at a single TE. Some of these studies utilized the linear TE-dependence of percent change of BOLD signal, which is a consequence of R2* signal decay (Menon et al., 1993). Barth et al. attempted to denoise data by incorporating the relationship between TE and signal change with a fuzzy cluster analysis (Barth et al., 2001). Peltier and Noll later used multi-echo fMRI to demonstrate that the percent signal changes of resting BOLD fluctuations demonstrate linear TE-dependence (Peltier and Noll, 2000, 2002).

One major application of multi-echo fMRI has been for BOLD contrast optimization by combining time courses of different TEs using a weighting scheme. Several weighting schemes have been proposed (Poser et al., 2006). Simple schemes use estimates of contrast-to-noise ratio as weights. More robust schemes use weights from contrast curves that are modeled for each voxel after estimating T2* from multi-echo fMRI data (Posse et al., 1999). Benefits of contrast optimization, as described by Posse et al., include reduction of susceptibility artifact and thermal noise. De-noising constitutes another application of multi-echo fMRI. Dual-echo fMRI de-noising approaches involve the regression of short TE, minimally BOLD weighted time series from longer TE, BOLD sensitive time series (Glover et al., 1996). Buur showed that both general least squares and ICA approaches weighted by TE-dependence factors could decouple the effects of collinear head motion from BOLD signal in multi-echo data (Buur et al., 2008, 2009). Most recently, the subtraction of early TE from late TE time series was shown to remove the effect of signal drift (Beissner et al., 2010). In summary, substantial information is uniquely available in multi-echo fMRI data that can be used to improve data quality.

Classification of ICA Components Based on TE-Dependence

In this study, resting state data were acquired with multi-echo fMRI to allow for differentiation of BOLD from non-BOLD signal components based simply on a goodness of fit to a R2* change model for multi-echo data. Components were first identified using decomposition of multi-echo data with ICA. Components were then analyzed to determine if signal changes were associated with R2* changes. Scores to summarize the overall component-level modulation as R2* change were computed. Sorting components according to these scores identified BOLD fluctuations without anatomical templates or time course priors and identified artifact fluctuations without time course models for motion or physiology. Removing non-BOLD fluctuations from time series data enabled data de-noising, which significantly improved seed-based functional connectivity measures, especially between subcortical and cortical regions.

Theory

Assuming mono-exponential decay, the signal across multiple echo times, TEn, where n is the echo number, varies as a function of initial signal intensity when the TE = 0 (S0) and relaxation-rate (R2*) according to Equation 1,

S(TEn) = S0exp(−R2 ∗ TEn),

(1)

where R2* is the inverse of relaxation time or 1/T2*. S0 can be modulated by changes in T1 (longitudinal relaxation rate), inflow, and motion. R2* varies as a function of magnetic field homogeneity, and specifically is modulated by changes in microscopic susceptibility due to changes in blood oxygenation. Figure 1 shows how changes in S0 (left column) and R2* (right column) affect the signal decay (top row) as a function of TE as well as the signal difference (middle row) and percent signal change (bottom row). Note that these difference and percent change curves can be used to differentiate whether or not the source of the signal change is BOLD (i.e. R2* change) or non-BOLD (i.e. S0 change) Expanding Equation 1 with a first order approximation for a small change in R2* gives:

Δ S/S = Δ S0/S0 − ΔR2 ∗ TE

(2)

An external file that holds a picture, illustration, etc. Object name is nihms363220f1.jpg

Shown are three echo simulations of BOLD (R2* change) and non-BOLD (S0 change) signals as a function of echo time (TE). The left column shows how the signal evolves for non-BOLD effects and the right column shows how the signal evolves for BOLD effects. The top row shows the signal during state x (no activation) and state y (activation). This top row demonstrates how the decay curves between rest and activation change in a different manner depending on if there is a change in (a) S0 or (b) R2*. The middle row shows the difference (y − x) signal for (c) change in S0, and (d) change in R2*. The bottom row shows the percent signal change (y−x)/0.5(x+y) for (e) change in S0, and (f) change in R2*.

Equation 2 shows that the linear relationship in Figure 1f is solely a function of TE and has slope equal to ΔR2*. The equation also shows that for mono-exponential decay and small changes in ΔR2*, the signal changes from modulations in S0 and R2* are linearly separable.

Multi-echo measurements can be acquired in fMRI using a single-shot, multiple gradient-echo EPI (echo planar imaging) sequence. For each time point in the fMRI time series, images are acquired at two or more different echo times (TEs). Figure 2a shows images from acquisition at 3 TEs. Figure 2b shows the three time series corresponding to each TE. From these three time series, the mean R2* and S0 are estimated from a fit to Equation 1.

An external file that holds a picture, illustration, etc. Object name is nihms363220f2.jpg

(a) Multi-echo EPI images acquired at TE values of 15ms, 39ms, and 63ms. Image intensity decreases exponentially with TE. (b). Left: multi-echo EPI time courses from a voxel in visual cortex (center of yellow box in (a)) during periodic visual stimulation plotted as percent signal change. Right: percent signal change amplitude as a function of TE (black), with linear fit, i.e. change in R2* (gold). The fit is significant (p<0.01), with ΔT2*=0.3ms. (c) Left: multi-echo EPI time courses from a single precuneus voxel (center of white box in (a)) during rest, plotted as percent signal change. TE is indicated by the color. Right: percent signal change amplitude as a function of TE (black), with linear fit, i.e. change in T2* (gold). The fit is significant (p<0.01), with dT2*=0.3ms.

The changes from mean R2* and S0 that underlie a signal fluctuation can also be estimated. To estimate changes in mean R2* and S0, a reference function, such as for a task design or resting fluctuation, is first regressed to the signal time course of each TE. TE-specific signal changes can be fit to the TE-dependence model in Equation 2, which estimates both ΔR2* and ΔS0 with one least squares fit. Fitting both ΔR2* and ΔS0, rather than one at a time, is unstable (Gowland and Bowtell, 2007). Therefore we separate Equation 2 into two sub-models, one for estimation of ΔR2* and one for ΔS0:

ΔS/S = Δ S0/S0 and ΔS/S = −ΔR2 ∗ TE.

(3)

Fitting signal changes to these sub-models has greater stability than fitting for both parameters simultaneously. This approach is critical for multi-echo acquisition when signals are acquired at a small number of TEs in the presence of the typical noise levels in fMRI. Some physiological or motion artifacts may produce coupled R2* and S0 changes, in which case the precise values of R2* and S0 signal change cannot be computed using this approach (Wu, 2005). However, for the purposes of classifying a fluctuation as mainly an R2* fluctuation and not an S0 fluctuation, the proposed approach is sufficient.

Goodness of fit statistics can be computed for the fits to these TE-dependence models. This computation can be performed voxel-wise with an F-test comparing the residual from the fit of a model to the residual of the null (zero) model (equal to the sum of the squares of signal changes). A separate F-value is computed for the ΔR2* sub-model and for the ΔS0 sub-model. F values are computed for 1 degree of freedom used and n-1 degrees of freedom remaining, where n is the number of echoes. The cumulative distribution function for F(1,n-1) can be used to compute p-values corresponding to F-values. For a given time series reference function, the maps of ΔR2* and ΔS0 can be thresholded according to these p-values. Furthermore, these F-values can be averaged, weighted by total signal power

where i is the TE index, n is the total number of echoes, and ΔSTE_i is the coefficient of the reference function and the time course at TEi. This produces two summary statistics, κ and ρ,

where v is the voxel index, m is the number of voxels in the brain. κ and ρ reflect the goodness of fit to ΔR2* and ΔS0 models respectively and convey a representative F value for the voxels with the largest signal changes. F-values are weighted by signal power so that κ and ρ are less representative of F-values for the small component signal changes, which are more affected by ICA estimation error. κ and ρ are used to rank how well components of linear models (here corresponding to ICA component time courses) agree with signal changes described by ΔR2* and ΔS0 signal models.

Methods

Subjects

Nine right-handed healthy volunteers participated in the study (7 male, 2 female). Informed consent was obtained under an approved National Institute of Mental Health protocol.

Data Acquisition

Imaging was performed on a General Electric (GE) 3 Tesla Signa HDx MRI scanner (Waukesha, WI). The scanner’s body coil was used for RF transmission, and an 8-channel receive-only head coil (GE, Waukesha, WI) was used for signal reception. High-order shimming was performed to minimize field inhomogeneity.

Anatomical images were acquired using a T1-weighted MPRAGE sequence (FOV 240mm, 224×224 in-plane resolution, TI 725 ms, SENSE (GE ASSET) acceleration factor 2). Functional images were acquired with a multi-echo EPI sequence (TR 2.5 s, flip angle 90, matrix size 64×64, in-plane resolution 3.75 mm, FOV 240 mm, 31 axial slices, slice thickness 4.2 mm with 0.3mm gap, acceleration factor 2). Three echoes were acquired with the shortest possible echo times, TE = 15ms, 39ms, and 63ms. The readout window width for each image was 24 ms, and receiver bandwidth was 125 KHz. Images were reconstructed off-line using a C implementation of SENSE (Pruessmann et al., 1999). A separate fast gradient echo scan (TR 150 ms, TE 2.1ms) with the same coverage as the multi-echo acquisition was used for SENSE calibration. The subjects were instructed to rest with eyes open and fixate on a cross-hair. Two resting functional runs of 148 images (time series duration = 6 min 10 sec) were acquired. Pulse and respiration data were acquired using scanner-integrated photoplethysmograph and respiratory bellows.

Data Pre-processing

Data pre-processing was performed using AFNI (Cox, 1996). Each of the steps mentioned below includes the corresponding AFNI function italicized in parentheses. Each functional run was pre-processed separately as follows. For the standard analysis pathway, RETROICOR corrections were applied first, based on the pulse and respiration data (3dretroicor) (Glover et al., 2000). For the ME-ICA analysis pathway, the following preprocessing was performed on unprocessed time series data. The first four time points were discarded to allow for magnetization to reach steady state. Slice time correction was applied (3dTShift). Motion correction parameters were estimated for each time point by aligning the middle TE (39ms) images to corresponding first time point image using a rigid-body (6 parameter) alignment procedure (3dvolreg). The functional to structural co-registration parameters were estimated by registering the skull-stripped middle TE image from the first time point to the skull-stripped anatomical image using an affine (12 parameter) alignment procedure with the local Pearson correlation cost-functional (3dSkullStip, 3dAllineate) (Saad et al., 2009). Motion correction and anatomical co-registration parameters were then applied in one step (3dAllineate). A brain mask was computed from the mean image of the shortest TE (15ms) time series (3dskullstrip) and applied to all images. Each image was spatially smoothed with a 5 mm FWHM Gaussian kernel within the functional mask (3dBlurInMask). Finally, all voxel time series were high-pass filtered for frequencies above 0.02 Hz.

The processing sequence branches at this point. The above-mentioned pre-processed data was used in the following two ways. First, as is typical for multi-echo fMRI data, the three echoes were combined to form a single time series. Regression analysis was performed on this data. Second, to de-noise this data, nuisance regressors were obtained from our novel multi-echo (ME)-ICA component sorting method described below.

Multi-echo Combination

Following the above preprocessing steps, we combined multi-echo time courses with a T2* weighting scheme (Posse et al., 1999) to produce a single time series data set. This weighting scheme enabled the application of the same weights to raw data and data denoised with different strategies. For each voxel, the mean of each of the three time courses was computed and used to estimate the overall baseline S0 and baseline T2* by fitting to Equation 1 with log-linear regression. The time courses were then optimally combined by weighted summation by a factor, w, described in Equation 6.

w(T2∗)n=TEn·exp(−TEn/T2(fit)∗)∑nTEn·exp(−TEn/T2(fit)∗).

(6)

On this data set, functional connectivity analysis using a seed voxel approach was performed using either standard de-noising or ME-ICA de-noising, as described below.

Nuisance Regressor Selection with ME-ICA

To find time series that are common across both TE and spatial location, TE was treated as a fourth spatial variable for spatial ICA. The data were decomposed as described in Equation 7:

D(Nt, Nx x Ny x [Nz x Ne]) = M(Nt, Nc) x C(Nx x Ny x [Nz x Ne])

(7)

Where D are the time series data, Nx, Ny, and Nz are the coordinates of voxels within the functional brain mask, Ne is the number of echoes, Nt is the number of time points, M is the mixing matrix (i.e. ICA component time courses), Nc is the number of components, and C are the component maps. To implement this spatial decomposition using MELODIC ICA (Beckmann and Smith, 2004b), the multi-echo data were concatenated over space. For each time point, the three volumes corresponding to each TE were combined into a single volume of size Nx, Ny, [Nz × Ne] (3dZcat). Each time course was de-meaned and variance-normalized. Dimensionality was then automatically estimated using probabilistic PCA, followed by dimensionality reduction to the estimated number of components (Beckmann and Smith, 2004a). FastICA (Hyvarinen, 1999) was applied to dimensionally reduced data to produce the mixing matrix, M.

For each ICA component, the TE-dependence of the component time series was mapped to localize where the component represented an S0 or an R2* modulation (as described in the Theory section). For each voxel, the TE-specific signal changes for the ICA component and TE-specific signal means were fit to compute ΔR2* or ΔS0:

ΔS = Δ S0/S0 ∗ S and ΔS = −ΔR2 ∗ TE ∗ S

(8)

using weighted ordinary least squares. Equation 8 is the formulation of Equation 2 under the assumption that thermal noise variance is independent of TE. For each sub-model, fit parameters, goodness of fit statistics F{ΔR2*} or F{ΔS0}, and p-values were mapped. Alpha probability simulations were computed for p(F)<0.05, and used to threshold the maps at p<0.05, cluster corrected for α<0.01. To summarize the overall ΔR2* and ΔS0 effects of the component, the averages of F weighted by total power of signal changes (sum of squares) were computed. κ represents the weighted average of F{ΔR2*} and ρ represents F{ΔS0}. High κ indicated strong ΔR2*-like character (BOLD-like), and high ρ indicated strong ΔS0-like character (non BOLD-like).

The ICA components were rank-ordered based on their κ and ρ scores. These two rank orderings (κ-spectrum and ρ-spectrum) were used to differentiate BOLD components from non-BOLD components. Both κ and ρ spectra were found to be L-curves with well-defined elbows distinguishing high score and low score regimes. This inherent separation was used to identify BOLD components in an automated procedure. First, the elbows of κ and ρ spectra were identified. The spectra were scanned from right to left to identify an abruptly increased score following a series of similarly valued low scores. The κ and ρ scores marking abrupt changes were used as thresholds. Those components with κ greater than the κ threshold and ρ less than the ρ threshold were considered BOLD components. All other components were considered non-BOLD components. These were used as noise regressors in time course de-noising.

De-noising and Functional Connectivity Analysis

De-noising by removing the nuisance regressors identified by ME-ICA was compared to de-noising by subtracting the fits of standard noise regressors and band pass filtering. Standard noise regressors included polynomial drifts, motion parameters, RETROICOR (Glover et al., 2000) and RVT (respiration variation over time) (Birn et al., 2006). These were fit to optimally combined time courses using multiple regression. Residuals were band pass filtered to a frequency range of 0.02 Hz to 0.1 Hz to produce data for functional connectivity analysis. The number of degrees of freedom used for de-noising was calculated as the number of noise regressors, excluding drift terms, plus the number of Fourier terms removed in the band pass filtering.

For seed-based connectivity analysis with standard de-noising (i.e. motion parameters, RETROICOR, RVT, and band-pass filtering), a seed time course was first extracted for a region of interest from the de-noised and band passed data. The seed time course was then regressed to all other de-noised and filtered time courses in a multiple regression model. To properly account for the degrees of freedom used in de-noising and band-pass filtering, a zero-column for each degree of freedom was included in the regression model.

For ME-ICA de-noising, nuisance regressors (low κ, high ρ) were removed from optimally combined data by first fitting the multiple regression model including all component time courses (the mixing matrix) and polynomial drifts, then subtracting the partial fit of the nuisance regressors and drifts to produce ME-ICA de-noised time courses. This is equivalent to the component removal procedure implemented in the FSL function fsl_regfilt. No other filtering was applied. The number of degrees of freedom used in de-noising was the number of nuisance regressors removed plus the number of drifts removed.

For seed-based connectivity analysis with de-noised ME-ICA data, a seed time course was extracted for a region of interest from ME-ICA de-noised data. This was regressed to all other de-noised, optimally combined time courses in a multiple regression model that included a zero-column for each degree of freedom used. R2 and T statistics were computed for all fits. T values accounted for degrees of freedom used. Group analysis of seed-based connectivity was done by transforming R values to Z values using the Fischer Z-transform, warping Z maps to Talairach space (auto_tlrc), and performing a one-sample T-test (3dttest).

Results

TE-dependence of BOLD-like and non BOLD-like components

The TE-dependence mapping is demonstrated for a BOLD-like and non BOLD-like related component in Figure 3. The BOLD-like component has a low-frequency time course and high percent signal change localized to the middle frontal and inferior parietal regions. There is clear correspondence between the localization of high percent signal changes, high ΔR2*, and strong goodness of fit to the ΔR2* model, F{ΔR2*}. Thresholding ΔR2* by F{ΔR2*} for p<0.05 and alpha<0.01 cleanly localizes BOLD-like signal fluctuations. In contrast, there is no correspondence between the localization of percent signal change and either ΔS0 or goodness of fit to the ΔS0 model, F{ΔS0}. Thresholding S0 by F{ΔS0} for p<0.05 and alpha<0.01 results in an empty map. κ for the BOLD-like component was 184, and ρ was 15.

An external file that holds a picture, illustration, etc. Object name is nihms363220f3.jpg

TE-dependence maps of ICA components from ME-ICA. (a) A BOLD-like and (b) a non BOLD-like component. For each component (a and b), the left panel shows percent signal change maps for three TEs 15ms, 39ms, 63ms (above), and the component time course (below). The right panel shows results of fitting to the ΔT2* change model (above) and the S0 change model (below). Goodness of fit maps, F{ΔT2*} and F{ΔS0}, are used to threshold parameter maps α<0.01 (p<0.05). (a) BOLD-like component: High percent signal change in gray matter scales linearly with TE. The component time course exhibits low frequency fluctuations. (b) non BOLD-like component: High percent signal change at the edge of brain is constant with TE. The component time course exhibits high frequency fluctuations.

The non BOLD-like component has a high frequency time course with localization along the brain edges. TE-dependence mapping of the artifact contrasts with TE-dependence mapping of the functional network. There is little correspondence between the localization of large percent signal changes and high ΔR2*, and no correspondence with strong goodness of fit to the ΔR2* model. Thresholding ΔR2* by F{ΔR2*} results in an empty map. There is good correspondence between localization of high percent signal change, percent S0 change, and goodness of fit to the S0 model. Thresholding ΔS0 by F{ΔS0} clearly localizes the artifact to brain edges. κ for the artifact component was 22, and ρ was 90.

Ranking ME-ICA components by κ

Figure 4a shows κ score vs the rank according to variance explained (i.e. ICA rank). This shows that the correspondence between these two measures is low. Figure 4b shows κ score vs. the rank according to κ. The rank ordering of κ scores, (the κ spectrum), show a clear L-curve behavior with distinct high κ and low κ regimes. Figure 4c shows that this distinct behavior is highly reproducible across subjects. Low κ values had mean of 21.9±0.5. High κ scores had a mean of 91.5±2.9. Figure 4d shows that the maps corresponding to the high κ scores appear to be similar to previously identified resting state network components (Damoiseaux, 2006b). The thresholded ΔR2* maps are shown for the top 12 components as ordered by κ scores. Representative high κ components include the default mode (IC 23), the sensory network (IC 13), and the motor network (IC 35). BOLD-like component time courses are not necessarily low frequency (< 0.1 Hz). For example, the time course for IC 30 has a higher κ score than the time course for IC 13.

An external file that holds a picture, illustration, etc. Object name is nihms363220f4.jpg

For a representative subject, κ score vs (a) ICA rank (variance explained), and (b) rank by κ (κ spectrum). The κ spectrum, is an L-curve with two distinct regimes: high κ (κ >20) and low κ (κ <20), with low κ components on a linear tail. (c) κ spectra for 8 subjects. (d) First 12 ME-ICA components ranked by κ for a representative subject. Each panel shows the time course and thresholded ΔR2* map. Components are annotated with κ-score, ρ-score, and ICA component number. All high κ components are clearly functional networks.

Ranking ME-ICA components by ρ

Figure 5a shows ρ score vs the rank according to variance explained. This shows that the correspondence between these two measures higher than that of κ with variance explained. Figure 5b shows ρ vs. the rank according to ρ. The rank-ordering of ρ scores, (the ρ spectrum), show a clear L-curve behavior with high ρ and low ρ regimes. Figure 5c shows that this distinct behavior is highly reproducible across subjects. Low ρ scores had mean of 24.3±0.8. High ρ scores had a mean of 53±3.1. Figure 5d shows thresholded ΔS0 maps for the top 8 components as ordered by ρ scores. Some of the representative high ρ components have high ΔS0 at brain edges (IC 3, 18, 4) and appear to be motion related. Some non BOLD-like component time courses have low frequency (> 0.1 Hz) contributions.

An external file that holds a picture, illustration, etc. Object name is nihms363220f5.jpg

For a representative subject, ρ score vs (a) ICA rank (variance explained), and (b) ρ rank (ρ spectrum). The ρ spectrum, like the k-spectrum, is an L-curve with two distinct regimes: high ρ (appx. ρ >20) and a linear tail with low ρ (appx. ρ <20). (c) ρ spectra for 8 subjects. (d) First 8 ME-ICA components ranked by ρ for a representative subject. Each panel shows the time course and thresholded %ΔS0 map. Components are annotated with κ-score, ρscore, and ICA component number. All high ρ components are clearly artifacts.

Components at the Elbow

Most components were unambiguously classified into high and low regimes for κ and ρ using the elbows of the corresponding spectra. Inevitably, a number of κ values were near the cutoff threshold (the elbow of the ranking curve). Figure 6 shows the TE-dependence maps for the ΔR2* and ΔS0 models for two of these components. IC 2 localizes high ΔR2* to the Circle of Willis and surrounding areas, and IC 17 localizes high ΔR2* to white matter and/or cerebral spinal fluid (CSF). The time courses of these components were low frequency. Comparing component time courses to RVT regressors showed strong correlation with two of the standard RVT regressors. This indicates that a component with a near-threshold κ score could reflect ΔR2* modulation from respiratory variation or related BOLD-like effects of no interest. The component with high ΔR2* near the Circle of Willis and surrounding area, IC 2, was selected as a regressor of no interest on the basis of having a high ρ score. The component with high ΔR2* in white matter/CSF was not selected as such due having a low ρ score.

An external file that holds a picture, illustration, etc. Object name is nihms363220f6.jpg

Components with κ scores near κ thresholds are correlated to low-frequency RVT time courses. Components are annotated with κ score, ρ score, and ICA component number. TE-dependence maps for ΔR2* and ΔS0 models show high ΔR2* localized to non-gray matter regions.

De-noising Resting Data in a Single Subject

Figure 7 shows the use of the ME-ICA nuisance regressors to de-noising the data. Figure 7a shows a comparison of time series from the corresponding single voxels in the center of the red boxes corresponding to the right insula, left hippocampus, and brain stem respectively. The top row of plots have only drifts removed. The second row shows the plots after drifts were removed, RETROICOR applied, and RVT and motion regressed out. The third row of plots shows the above time course with band-pass filtering applied. The fourth row of plots shows the time courses with drifts and ME-ICA derived nuisance time series (low κ and high ρ) regressed out. Time courses from these areas after removal of drifts show substantial high frequency noise and spiking. Standard de-noising by removing RETROICOR, RVT, and motion regressors reduces high frequency noise. These time courses are somewhat smoothed by band-pass filtering. ME-ICA de-noising reduces high frequency noise, spiking, as well as some low frequency fluctuations, without the use of physiological noise modeling or band-pass filtering.

An external file that holds a picture, illustration, etc. Object name is nihms363220f7.jpg

Signal from three regions of interest from a representative subject: the right insula, left hippocampus, and brainstem. (a) shows de-noising by removing: drifts only; drifts, physiology, and motion (the standard); drifts and low κ components (ME-ICA). (b) Seed based connectivity measured by R2 and T values, with baseline regression for: motion and physiology, then band pass filtering for 0.02–0.1 Hz; drifts and low κ components, without band pass filtering.

Figure 7b shows seed-based functional connectivity maps obtained using seed time courses from the de-noised data. R2 correlation maps show that ME-ICA de-noising, without band pass filtering, reveals greater functional connectivity to gray matter clusters than de-noising with standard noise regressors and band pass filtering. Axial views of R2 maps for insula and hippocampus connectivity show that the de-noising methods produce similar connectivity patterns proximal to the seed, but ME-ICA de-noising exposes greater long distance correlation. With ME-ICA de-noising, the insula shows greater correlation to premotor and cingulate regions, hippocampus shows greater correlation to premotor and sensory regions, and brainstem shows greater correlation to frontal and parietal regions. T-maps show that T-statistics are much higher for correlation with ME-ICA de-noising than for correlation with standard de-noising and band pass filtering.

Application to Group Level Correlation Maps

Group-level connectivity was evaluated using one-sample T-tests of the individual-level correlation maps from standard and ME-ICA based de-noising. Unthresholded group T-maps for hippocampus and brainstem connectivity are shown in Figure 8 for ME-ICA and standard de-noising. The group T-maps based on low κ de-noising showed much higher T-statistics for connected regions than the group T-maps based on standard de-noising. This indicated that (Z-transformed) correlation coefficients based on ME-ICA were more consistent across subjects than Z-transformed correlation coefficients based on standard de-noising. Comparing Figures 7 and ​8 show that for maps based on ME-ICA de-noising, the regions of higher group T-value correspond to the regions of higher individual level T-values from regression analysis (white arrows). Figure 9 shows thresholded axial views of the connectivity maps in Figure 8. The brainstem shows clear connectivity to putamen and caudate nuclei in the basal ganglia and to bilateral premotor and parietal cortical areas. The hippocampus shows clear connectivity to bilateral sensory (dorsal and ventral postcentral gyrus), temporal, and premotor cortical areas. Thresholding group-level connectivity maps that were based on standard de-noising produces empty maps at FDR corrected q<10−4, which is the significance level used in Figure 9.

An external file that holds a picture, illustration, etc. Object name is nihms363220f8.jpg

Group T-maps of subcortical connectivity with hippocampus and brainstem computed from individual maps with baseline regression for: motion and physiology, then band pass filtering; drifts and low κ components, without band pass filtering.

An external file that holds a picture, illustration, etc. Object name is nihms363220f9.jpg

Group maps for subcortical connectivity with hippocampus and brainstem after removal of drifts and low κ and/or high ρ components. Overlay is map of mean Z-value, thresholded by T-value corresponding to FDR corrected p<10−4. Underlay is template brain in Talairach space.

Discussion

The differentiation of BOLD and non-BOLD signal is a fundamental problem in resting state fMRI analysis. In the past, approaches to this problem have included regression of temporal models derived from motion and physiologic signals. Here, an approach is introduced that involves acquiring resting state data with multi-echo EPI, identifying BOLD-like (high κ, low ρ) non-BOLD-like (low κ, high ρ) components directly from the data, and using these non BOLD-like components to obtain nuisance regressors. This approach to selecting components and de-noising does not require external physiologic measures, temporal noise models, or anatomical templates, and is fully automated. It is based instead on ICA and the principle that the BOLD signals of resting neural activity are characteristically TE-dependent.

The current study benefits from previous research on TE-dependence and multi-echo fMRI. The TE-dependence of activation-induced signal changes has been demonstrated multiple times since 1992 (Ogawa et al., 1992; Bandettini et al., 1994; Menon et al., 1995). Speck and Henning mapped T2* and S0 for activation corresponding to visual stimulation (Speck and Hennig, 1998). De-noising with dual-echo acquisition has been demonstrated in several instances (Glover et al., 1996; Buur et al., 2008; Beissner et al., 2010, 2011). The TE-dependence of low-frequency resting state fluctuations was demonstrated with 4-echo acquisition (Peltier and Noll, 2002). Poser et al. introduced accelerated acquisition of multi-echo EPI using parallel imaging (Poser, 2006), which enabled the acquisition of whole brain multi-echo fMRI within a conventional TR.

Wu et al. had demonstrated that essentially all fMRI artifacts, regardless of their origins in motion, physiology, or other sources could be expressed in terms of R2*, S0, and combinations thereof (Wu and Li, 2005). The current method extends that framework to the classification of ICA components as R2*, S0 or coupled components. It also differs from previous approaches such as the estimation of ΔR2* and/or ΔS0 at each time point to produce parameter time courses (Speck and Hennig, 1998). Simulations demonstrated that point wise fits produce ΔR2* time courses with higher levels of noise than signal time courses of a single echo (Gowland and Bowtell, 2007). In our approach, regression weights were fit to ΔR2* and ΔS0 models, which is similar to the approach of Peltier and Noll (Peltier and Noll, 2002). The robustness of the fits in their study and in ours is attributed to the averaging effect of using a fluctuation over time to estimate ΔR2*. Our approach is also similar to that of Buur et al. (Buur et al., 2008), in which task-related BOLD signal was separated from co-linear head motion artifact by un-mixing using ICA with TE-dependence weights. Buur et al. concluded that combining ICA and linear TE-dependence analysis is more robust than other methods of extracting BOLD fluctuations from multi-echo data. Our results show that TE-dependence separates BOLD fluctuations from many kinds of non-BOLD artifacts, including both motion and physiology. Furthermore, the results observed after ME-ICA de-noising at the individual level are reflected at the group level.

fMRI artifacts arising via non-R2* mechanisms include S0 fluctuations as well as coupled R2* and S0 fluctuations such as RVT-like components. Non-R2* components were clearly distinguished as low-κ components, and specifically S0 fluctuations were further characterized with high ρ values. The intermediate κ and ρ values for RVT components were due to signal changes of these components fitting both R2* and S0 models. Wu et. al had previously suggested that some sources of fMRI artifact would produce coupled R2* and S0 changes (Wu and Li, 2005). The results showed that components of coupled origin could be identified by their rankings on κ and ρ-spectra.

Prior studies have denoised task-based fMRI by removing ICA components that have poor correlation to task reference functions (Thomas et al., 2002). Other studies have proposed de-noising resting state fMRI by removing those ICA components that do not correspond to atlases for “functional networks” (De Martino et al., 2007; Perlbarg et al., 2007). This is a questionable procedure because removing components that are neither networks nor obvious artifacts (McKeown et al., 1998; Beckmann and Smith, 2004b) depends on a circular argument. We show for the first time using ICA of multi-echo fMRI that network-like ICA components fit well to a BOLD TE-dependence model, which would be expected, but also that all other components fall into a well-defined non-BOLD regime, which is a novel and important finding. This separation was remarkably stable across subjects, enabling fully automated identification of network-like components and robust time course de-noising without spatial templates or time course models for noise.

Using data that was de-noised with ME-ICA nuisance regressors in seed-based correlation analysis resulted in significantly improved correlation maps in individual subjects and in groups. Studying functional connectivity of subcortical regions is challenging due to low functional contrast-to-noise due to CSF and blood flow pulsatility and distance from receiver elements. Where standard de-noising showed no clear correlation patterns for the hippocampal and brain stem seeds, ME-ICA de-noising revealed robust correlation patterns. The brain stem seed was localized to the anterior pons that contains corticospinal (pyramindal) tracts connecting to premotor, parietal, and motor regions (Kiernan, 2009). This pattern of anatomical connectivity agrees well with the pattern of functional connectivity exposed after ME-ICA de-noising. The hippocampus seed was localized to the head of the right hippocampus that has anatomical connectivity to sensory regions via temporal and entorhinal cortices (Kiernan, 2009). The pattern of functional connectivity exposed after ME-ICA denoising agreed with this pattern of anatomical connectivity. These enhancements can be attributed to the robustness of using information across both space and TE to extract components and then identifying component origin by evaluating goodness of fit to the BOLD TE-dependence model.

Overall, analysis of ΔR2* and ΔS0 for ME-ICA components was shown to be a powerful approach to differentiating BOLD and non-BOLD signal in resting state data for both the seed-based and ICA approaches to connectivity analysis. In particular, the proposed method shows considerable promise in removing the pulsatile artifactual signal in subcortical regions such that subcortical-cortical connectivity can be studied more effectively. The benefits also extend beyond the application to resting state data. The improved contrast to noise will have direct impact on the repeatability and quality of activation-related fMRI studies. This method will be particularly beneficial for clinical fMRI of patients who exhibit a high amount of movement. It also holds potential for reducing the effects of stimulus correlated motion as in studies of overt speech production.

In the present study, an 8 channel head coil, 3T scanner, and accelerated imaging (SENSE factor 2) was used to acquire whole brain images at 3 TEs with a 3.75 mm × 3.75 mm × 4.5 mm resolution and a 2.5s TR. Recent advances in coil technology (Wiggins et al., 2006) and multi-slice excitation (Feinberg et al., 2010; Moeller et al., 2010; Setsompop et al., 2011) should allow for substantial increases in resolution and efficiency. These advances will allow for acceleration factors of 3 in plane and 3 in the slice direction, enabling the acquisition of whole brain images at 3 TEs with a 2mm × 2mm × 2mm resolution and a 2s TR. Acquiring more than 3 TEs for the purposes of the proposed analysis is a subject of further study. Given the good differentiation of components at the current resolution with 3 TEs, the use of more TEs could be redundant. For lower SNR regimes that are associated with higher resolution acquisition, acquiring more TEs could be beneficial. Given a particular TR, however, increasing resolution competes with acquiring more TEs, so this trade-off defines a limitation of the proposed approach.

ME-ICA allows for robust assessment of resting state correlation, producing maps from individual subjects where multi-subject averaging was previously required. This increase in functional contrast to noise in time series will likely lead to more robust clustering and segmentation of individual subject resting state data (Cohen et al., 2008) and substantially improve the power of multi-subject studies (Biswal et al., 2010).

​

Highlights

Acknowledgments

This work was funded by the National Institute of Mental Health Intramural Research Program. Prantik Kundu is supported by the NIH-Cambridge Scholars program.

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