Global and System-Specific Resting-State fMRI Fluctuations Are Uncorrelated: Principal Component Analysis Reveals Anti-Correlated Networks

Abstract

The influence of the global average signal (GAS) on functional-magnetic resonance imaging (fMRI)–based resting-state functional connectivity is a matter of ongoing debate. The global average fluctuations increase the correlation between functional systems beyond the correlation that reflects their specific functional connectivity. Hence, removal of the GAS is a common practice for facilitating the observation of network-specific functional connectivity. This strategy relies on the implicit assumption of a linear-additive model according to which global fluctuations, irrespective of their origin, and network-specific fluctuations are super-positioned. However, removal of the GAS introduces spurious negative correlations between functional systems, bringing into question the validity of previous findings of negative correlations between fluctuations in the default-mode and the task-positive networks. Here we present an alternative method for estimating global fluctuations, immune to the complications associated with the GAS. Principal components analysis was applied to resting-state fMRI time-series. A global-signal effect estimator was defined as the principal component (PC) that correlated best with the GAS. The mean correlation coefficient between our proposed PC-based global effect estimator and the GAS was 0.97±0.05, demonstrating that our estimator successfully approximated the GAS. In 66 out of 68 runs, the PC that showed the highest correlation with the GAS was the first PC. Since PCs are orthogonal, our method provides an estimator of the global fluctuations, which is uncorrelated to the remaining, network-specific fluctuations. Moreover, unlike the regression of the GAS, the regression of the PC-based global effect estimator does not introduce spurious anti-correlations beyond the decrease in seed-based correlation values allowed by the assumed additive model. After regressing this PC-based estimator out of the original time-series, we observed robust anti-correlations between resting-state fluctuations in the default-mode and the task-positive networks. We conclude that resting-state global fluctuations and network-specific fluctuations are uncorrelated, supporting a Resting-State Linear-Additive Model. In addition, we conclude that the network-specific resting-state fluctuations of the default-mode and task-positive networks show artifact-free anti-correlations.

Introduction

S ince the pioneering work of Biswal and associates (1995), the study of spontaneous functional magnetic resonance imaging (fMRI) activity has emerged as an important field of research. The correlation observed in slow spontaneous blood oxygenation level dependent (BOLD) fluctuations (<0.1 Hz) between the left and right sensory-motor cortex (SMC) as well as between the left SMC and medial motor areas have been replicated in numerous studies. Moreover, it has been extended to observations of functional connectivity within several other resting-state networks (Biswal et al., 1995; Greicius et al., 2003; Fox et al., 2005; Fox and Raichle, 2007; Lowe et al., 1998). Here we investigated the relationship between resting-state fluctuations that are global, common to all resting-state networks, and network-specific fluctuations.

fMRI-measured resting-state fluctuations correlate with and are therefore likely to include contributions from slow, local fluctuations in neuronal activity (Shmuel and Leopold, 2008), and fluctuations in systemic physiological parameters (Beall and Lowe, 2007; Bianciardi et al., 2009; Birn et al., 2006; Lund et al., 2006; Wise et al., 2004). fMRI-measured resting-state fluctuations include a component that is common to the majority of networks, as well as network-specific components (Fox et al., 2009). The global component includes contributions from neurophysiological activity (Schölvinck et al., 2010) and systemic physiological changes (Beall and Lowe, 2007; Birn et al., 2006; Lund et al., 2006). Irrespective of its origin, the global component might obscure resting-state fluctuations that are more network specific. Indeed, removal of the global effect in fMRI resting-state studies has been shown to be useful in facilitating the detection of interactions that are known to exist at the neurophysiological level in animal models (Fox et al., 2009).

Previous PET and fMRI task/response-based studies characterized the global average signal (GAS), defined as an average of the imaging signal over the entire brain (Aguirre et al., 1998; Friston et al., 1990; Zarahn et al., 1997). Other response-based studies developed suitable strategies for its removal (Andersson et al., 2001; Desjardins et al., 2001; Gavrilescu et al., 2002; Macey et al., 2004). In task/response-based studies, two main modeling approaches have been considered to account for the effect of the GAS: the additive model and the proportional scaling model (Desjardins et al., 2001; Gavrilescu et al., 2002). The former assumes that a global signal is added on top of the stimulus-dependent brain activity, whereas the latter treats the GAS as a gain effect that influences all BOLD signals equally.

To account for global effects in the analysis of resting-state functional connectivity, the global signals as well as the network-specific fluctuations need to be correctly estimated. In resting-state fMRI studies, the common component has previously been estimated by the GAS, defined as the average of the BOLD signal over the entire brain (Fox et al., 2005). Regressing out the GAS and model-based removal of physiological artifacts enabled the detection of a negative correlation between resting-state BOLD signal fluctuations of the default-mode and the task-positive networks (Chang and Glover, 2009; Fox et al., 2005, 2009; Fox and Raichle, 2007; Fransson, 2005). Removal of the GAS may be inadvisable, however (Murphy et al., 2009; Weissenbacher et al., 2009), as it is prone to producing spurious correlation artifacts, introducing functional connectivities that do not exist in reality (Murphy et al., 2009). This undesired feature resulting from the regression of the GAS called into question the interpretation of negatively correlated resting-state networks (Murphy et al., 2009).

Here we aimed to find an estimator of the global fluctuations in resting-state fMRI studies. Specifically, we searched for an estimator immune to the complications associated with the currently used estimator, the GAS. We hypothesized that the global effect is uncorrelated with the network-specific fluctuations. To test this hypothesis, we applied principal component analysis (PCA) to spatio-temporal fMRI time-series of resting-state activity. We show that there exists a single temporal component (eigenvariate) that closely resembles the GAS and accounts for a large proportion of the variance. By means of numerical simulations, we demonstrate that this property is not mandated by the algebraic manipulations involved in PCA, but is intrinsic to the complex spatio-temporal structure of the resting-state fMRI data. By definition, this eigenvariate is uncorrelated to the remaining principal components (PCs) that represent the more network-specific resting-state fluctuations. We conclude that resting-state global and network-specific fluctuations are uncorrelated. We propose a novel estimator of the global effect in resting-state fMRI, defined as the PC that correlates best with the GAS. Irrespective of whether the fMRI-measured resting-state global effect does or does not reflect neurophysiological activity, the removal of the PC-based global effect estimator facilitates the observation of network-specific functional connectivity. By applying an artifact-free removal of our proposed PC-based global effect estimator, we show that resting-state fluctuations of the default-mode and the task-positive networks are anti-correlated. Preliminary results were previously reported in an abstract form (Carbonell et al., 2010).

Methods

Subjects

Data from human subjects were downloaded from the BS002 database (www.brainscape.org; Fox and Raichle, 2007; Fox et al., 2009). BOLD-sensitized fMRI data (4×4×4 mm voxels, TE 25 msec, TR 2.16 sec) were acquired from 17 normal right-handed young adults using a 3 Tesla Siemens Allegra MR scanner. All subjects completed four 7-min resting-state runs, during which 194 volumes were acquired. Subjects were instructed to fixate their eyes on a cross-hair, remain still, and not fall asleep. Structural data (for atlas transformation) included a high-resolution (1×1×1.25 mm) sagittal, T1-weighted MP-RAGE (TR 2.1 sec, TE 3.93 msec, flip angle 7°) and a T2-weighted fast spin echo scan.

Preprocessing of fMRI data

The fMRI data were preprocessed using the standard stereotaxic fMRI preprocessing pipeline implemented in the neuroimaging analysis kit (NIAK^*). The first three volumes of each run were discarded to allow the magnetization to reach equilibrium. Each dataset was corrected for inter-slice differences in acquisition time, rigid body motion, ultra-slow time drifts (high-pass filter with a 0.01-Hz cut-off), high temporal frequencies (low-pass filtering with a 0.1-Hz cut-off), and physiological noise (CORSICA; Perlbarg et al., 2007). For each subject, the mean motion-corrected volume of all runs was co-registered with an individual T1 scan using Minctracc (Collins et al., 1994), which was itself nonlinearly transformed to the Montreal Neurological Institute (MNI) nonlinear template using the CIVET pipeline (Zijdenbos et al., 2002). The functional volumes were re-sampled in the MNI space at a 2-mm isotropic resolution and spatially smoothed with a 6-mm isotropic Gaussian kernel. Regression of the time-courses of the six parameters resulting from rigid body motion correction was applied. In addition, each voxel time course was centered on zero by subtracting its mean value (over time) and was normalized by dividing by its root sum of squares (unit variance normalization).

For each individual run of fMRI time series, voxels in the brain were defined by constructing a brain mask using the following procedure. After motion correction, the mean fMRI volume was blurred using an isotropic Gaussian kernel (Full Width Half Maximum [FWHM] of three times the voxel size). An intensity-based threshold was then applied to derive a binary mask of the brain. The threshold value was estimated from the histogram of the blurred volume using the Otsu's algorithm (Otsu, 1979) for finding an optimal cut-off point for separating the mixture into two independent distributions (namely, the intensities of the background and of the brain).

Definition of the PCA-based global effect estimator

Let the BOLD time series be represented as the V×T matrix X , where V and T denote the number of voxels and time-dependent volumes of data, respectively. All resting-state fMRI studies that regress out the GAS implicitly assume that the BOLD measured resting-state signals include a global effect component that needs to be removed by standard regression techniques. Here, we explicitly formulate this assumption by introducing a Resting-State Linear Additive Model (RSLAM), in which a global effect signal is added on top of system-specific fluctuations. Thus, at any voxel v, this RSLAM can be expressed as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*}{ \bf X} ( v ) = { \bf g} ( v ) + { \bf X}_{RS} ( v ), \end{align*} \end{document}

Where g is identical in all the voxels up to a multiplicative amplitude coefficient (hence the term “global signal”), and X _RS is a linear mixture of network-specific components, with nonzero amplitude only in portions of the brain. A common procedure is to estimate the global effect g at all time points by the GAS \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$${ \bf g}_{Av}$$ \end{document} over all intra-cerebral voxels: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$ { \bf g } _ { Av } = \frac { 1 } { V } \sum \nolimits_ { v = 1 } ^V { \bf X } ( v )$$ \end{document} . A major drawback of this procedure is that it forces the spatial average of the resting-state fluctuations X _RS to be zero, and therefore it might mandate spurious negative correlations in seed-based correlation analysis (Murphy et al., 2009).

Here we use PCA (Baumgartner et al., 2000; Friston et al., 1993) to decompose the data into a set of orthogonal components. We show that the GAS \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$${ \bf g}_{Av}$$ \end{document} can be approximately captured by a single PC. PCA decomposes X into a linear combination of orthogonal components. A relatively small number of these components account for most of the variance in the fMRI data. This decomposition is given by the expression \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*}{ \bf X} = l_1{ \bf e}_1{ \bf s}_1 {\prime} + \ldots + l_r { \bf e}_r { \bf s}_r {\prime} , \qquad \qquad &(1), \end{align*} \end{document}

where the orthogonal vectors \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$${ \bf e}_1 , \ldots , { \bf e}_r$$ \end{document} are termed eigenimages, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$${ \bf s}_1 , \ldots , { \bf s}_r$$ \end{document} are the corresponding eigenvariates (temporal components), and the number of components r≤T is determined by the rank of the matrix X . In addition, the values l₁ ≥…≥l_r enable the quantification of the proportion of the total variance (PTVar) accounted for by each PC. The proportion of the total variance explained by the i-th component is given by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$PTVar ( PC_i ) = l_i^2 / \sum \nolimits_{k = 1}^r l_k^2.$$ \end{document}

By taking average in (1), one can demonstrate that there exists a vector \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$${ \bf a} = ( { \bf a}_1 , \ldots , { \bf a}_r ) ( { \bf a}^ \prime { \bf a} = 1 )$$ \end{document} such that the normalized GAS ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$${ \bf g}_{Av} / \sqrt{{ \bf g} \prime_{Av} { \bf g}_{Av}}$$ \end{document} ) can be expressed as a weighted linear combination of the eigenvariates, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$${ \bf a}_1{ \bf s}_1 \prime + \ldots{ \bf a}_r { \bf s}_r \prime$$ \end{document} Thus, the proportion of the total variance explained by the GAS \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$${ \bf g}_{Av}$$ \end{document} is given by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$PTVar ( { \bf g}_{Av} ) = \sum \nolimits_{k = 1}^r { \bf a}_k^2l_k^2 / \sum \limits_{k = 1}^r l_k^2$$ \end{document} .

As pointed out by Andersson and associates (2001), an ideal global effect estimator should explain more variance shared by virtually all voxels than any other global effect estimator. According to the properties of the PCA decomposition, the maximal proportion of the total variance explained by a weighted linear combination of the eigenvariates is reached when a is given by the canonical vector a \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$= ( 1 , 0 , \ldots 0 )$$ \end{document}

Our analysis consists of taking those canonical vectors of the form \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$${ \bf a}_i = ( 0 , \ldots , 0 , 1 , .0 )$$ \end{document} (1 in position i and 0 elsewhere, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$i = 1 , \ldots , r$$ \end{document} ) and finding the index \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$\hat{{ \bf i}}$$ \end{document} that maximizes \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$r_i = corr ( { \bf s}_i , { \bf g}_{Av} )$$ \end{document} . Then, our PC-based estimator of the global effect g is based on the single PC with the index \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$\hat{{ \bf i}}$$ \end{document} and is given by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$ { \bf g_ { PC } } = l_ { \hat { i } } ( \frac { 1 } { V } \sum \nolimits _ { v = 1 } ^V e_ { \hat { i } } ( v ) ) s_ { \hat { i } } \prime$$ \end{document} . Note that this definition of g _PC includes spatial averaging that scales the single PC estimator to amplitude similar to that of the global average.

Our proposed estimator has a number of appealing properties. First, as will be shown in the Results section, this time course shows high correlation with the GAS. It is therefore a good approximation of the global effect. Second, the condition \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$\hat{i} = 1$$ \end{document} (as it holds in most typical cases) would imply that \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$PTVar({\bf g}_{PC}) \ge PTVar ({\bf g}_{Av})$$ \end{document} , which by definition makes our estimator a better estimator of the global effect compared to the GAS (because it explains more variance; Andersson et al., 2001). Third, the decomposition (1) allows us to measure the local (voxel-wise) contribution [ g _PC (v)] of every voxel to the PC-based global effect estimator: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$g_{PC} ( v ) = l_{ \hat{i}}e_{ \hat{i}} ( v ) s_{ \hat{i}} \prime$$ \end{document} . In addition, due to the orthogonality property among the PCs, we obtain a global effect estimator that is orthogonal to the less-global, more-network-specific resting-state fluctuations, which can be estimated as the sum of the remaining components. Thus, a relatively high correlation between our PC-based global effect estimator and the GAS will support the assumed additive model for resting-state BOLD fMRI signals (due to the orthogonality property of PCA).

Spatial distribution of global effect

The significance of the contributions of different parts of the brain to the two global effect estimators considered here was estimated within the framework of General Linear Model (GLM). Specifically, we employed the model \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$${ \bf X} ( v ) = \beta_g ( v ) { \bf g + e}$$ \end{document} where g denotes either the GAS g _Av or the PC-based global effect estimator g _PC, and e is a vector of independent and identically distributed Gaussian variables. To detect regions that contributed significantly to the global effect g , we tested the null hypothesis that the coefficient \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$\beta_g ( v )$$ \end{document} is zero. This yielded the t Statistical Parametric Map (SPM) \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$T = {\hat \beta}_g ( v ) / std ( {\hat \beta}_g ( v ) )$$ \end{document} , where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$std ( {\hat \beta}_g ( v ) )$$ \end{document} denotes the standard deviation of the estimated coefficient \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$${\hat \beta}_g ( v )$$ \end{document} . In addition, the GLM framework allowed us to test the null hypothesis of no differences between the contributions of different parts of the brain to the two global effect estimators (i.e., \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$\beta_{g_{PC}} ( v ) = \beta_{g_{Av}} ( v )$$ \end{document} ). This test, too, yields a t SPM for the contrast \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$\beta_{g_{PC}} ( v ) - \beta_{g_{Av}} ( v )$$ \end{document} , which is given by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$T = ( {\hat \beta}_{g_{PC}} ( v ) - {\hat \beta}_{g_{Av}} ( v ) ) / \sqrt{std ( {\hat \beta}_{g_{PC}} ( v ) ) ^2 + std {\hat \beta}_{g_{Av}} ( v ) ) ^2}$$ \end{document} .

In this article, we used the PCA implementation given by the fMRIStat matlab toolbox.^† The parameter estimation in the framework of the SPM was also implemented using the fMRIStat matlab toolbox,^† including the correction for multiple comparisons (Worsley et al., 1998) taken into account while determining the threshold level of the t SPM maps. Further statistical analysis included the combination of runs within subjects and group analysis of all subjects in the sample. For this, we applied hierarchical random effect analysis proposed by Worsley and associates (2002) and implemented in the fMRIStat matlab toolbox.

Seed-based correlations: regressing PC-based estimator does not cause spurious correlations

To compare the effect of removing each of the two global signal estimators (the GAS and the PC-based global effect estimator) on seed-based correlation analysis of resting-state BOLD fluctuations, we applied two 6-mm spherical seeds. The seeds were positioned within the posterior cingulate/precuneus cortex (PCC; Talairach coordinates [-2-36 37]) and left middle temporal (MT) cortex ([-47-69-3]) (these spatial locations were also used by Fox et al., 2009). Seed-based t SPM of single runs were computed within the GLM framework by considering the seed time course of interest ( f ) as a covariate in the model and the global effect g (either GAS or PC-based estimator) as a confounding variable, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} X ( v ) = \beta_g ( v ) g + \beta_f ( v ) f + e \tag{2} \end{align*} \end{document}

Under the RSLAM assumption, model (2) can be rewritten as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} c ( v ) g + { \bf X}_E ( v ) = \beta_g ( v ) g + \beta_f ( v ) ( cg + f_E ) + e , \tag{3} \end{align*} \end{document}

where X_E(v) and f_E denote the system-specific components of X(v) and f, respectively, and c(v), c are multiplicative constants. Within this framework, regressing out the global effect would enhance the network-specific seed-based correlations (specified by β_f ). However, a critical issue here is the estimation of the global effect g . As pointed out by Murphy and associates (2009), estimating g as the GAS might introduce artificial seed-based correlations.

For analyzing the causes for the modification of seed-based correlations following the regression of the GAS, we separated this regression to two stages that contribute distinct modifications: modeling and estimation. Regarding the modeling stage, the assumption of the RSLAM and the subsequent regression of the global effect reduce the seed-based correlation values compared to the correlation values that could be obtained under a more simplistic model where no global effect is considered. Therefore, we can consider this overall decrease in the seed-based correlation values as merely the result of the assumed RSLAM. In other words, the appearance of certain modified correlations after the regression of the global effect g is completely justified by the assumption of the RSLAM. Regarding the estimation stage, a second source of spurious correlations arises at the estimation level, when estimating the global effect as the GAS. As we will show below, the GAS causes spurious correlations beyond the global connectivity strength reduction justified by the RSLAM.

Regressing out the global effect g from model (3) is equivalent to multiplying each term of the model by the projection matrix Pg=I-g(g'g)-1g’, which yields \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} P_gc ( v ) g + P_g { \bf X}_E ( v ) = \beta_g ( v ) P_gg + \beta_f ( v ) ( P_gcg + P_gf_E ) + e , \tag{4} \end{align*} \end{document}

Since P_gg=0, we then obtain from (4) the following model \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} P_g{ \bf X}_E ( v ) = \beta _f ( v ) P_gf_E + e , \tag{5} \end{align*} \end{document}

which is equivalent to model (6): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*}{ \bf X}_E ( v ) = \beta_g ( v ) g + \beta_f ( v ) f_E + e. \tag{6} \end{align*} \end{document}

Notice that in model (6), the global effect g plays the role of a confounding variable. Thus, the main challenge here is to quantify the impact of the global effect on the seed-based correlations as compared to the situation where no global effect is considered in the model. In other words, we need to quantify the impact of g on the coefficient β_f(v) estimated from model (6) as compared to the coefficient estimated from the model \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*}{ \bf X}_E ( v ) = \beta_f ( v ) f_E + e. \tag{7} \end{align*} \end{document}

Based on (Frank, 2000), the ordinary linear least squares estimator of the coefficient β_f(v) in (6) is given by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} \begin{split}& \beta_f ( v ) = \\ & \frac{ ( { \bf X}_E ( v )^ {\prime} { \bf X}_E ( v ) ) ^{1 / 2}} {( f_E^ \prime f_E ) ^{1 / 2}} \frac{corr ({ \bf X}_E ( v ) , f_E ) - corr ( { \bf X}_E ( v ) , g ) corr ( f_E , g )} { 1 - corr ( f_E , g ) ^2}.\end{split} \tag{8} \end{align*} \end{document}

Notice that the estimated coefficient in model (7) can also be obtained from expression (8) by setting g= 0. Hence, the inclusion of a global effect g satisfying the condition corr( g,f_E )≠0 in model (6) would modify the value of the estimated coefficient as compared to the estimated coefficient resulting from model (7), which represents the estimation of the system-specific correlation. In contrast, a global effect g satisfying the condition corr( g,f_E )=0 would not modify the original value of the estimated coefficient.

When g=g_Av (GAS), the global effect g and the seed time course f_E are correlated (corr( g,f_E )≠0) because by definition g includes the term f=g+f_E . Thus, regressing out the GAS introduces artificial seed-based correlations [modifications of the estimated coefficient in (7)]. In contrast, by construction, in the case g=g _PC (PC-based estimator) we have corr( g,f_E )=0 and no modifications to the original estimator (artificial correlations) are introduced. Notice that, in contrast to the proof given above, Murphy and associates (2009) used the spatial sum of the beta coefficients β_f(v) as an indicator for the possible appearance of spurious correlations following the removal of the GAS.

Finally, to test the significance of seed-based correlations we examine the null hypothesis β_f(v)=0, which yields a t SPM given by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document}$$T = {\hat \beta}_f ( v ) / std ( {\hat \beta}_f ( v ) )$$ \end{document} In a similar way, we also test the null hypothesis of no differences between the connectivity maps corresponding to the removal of each global effect estimator.

Numerical simulations

Numerical simulations were carried out to test whether the hypothesized high correlation between the GAS and a single PC is mandated by the algebraic manipulations involved in PCA. Our simulated volumes consisted of smooth Gaussian random fields (zero mean and standard deviation 1) sampled over a grid of 128×128×64 voxels of dimension 1×1×1 mm each. Here, each Gaussian random field was simulated by convolving white noise sampled in the 128×128×64 grid with a 3-dimensional Gaussian kernel of spatial FWHM equal to 6 mm in each direction (this value was chosen to mimic the spatial characteristics of measured fMRI data following preprocessing).

Four different types of datasets were generated, each sample consisting of a time series composed of 256 volumes (TR=1 sec). The first dataset type corresponded to independent (temporally uncorrelated) volumes. In contrast, for the other three types of datasets, a temporal correlation structure was added to the time series of 256 independent Gaussian random fields. In each case, an order 1 autoregressive model (AR1) was used to simulate temporally correlated volumes. Therefore, given an initial volume X₀ (Gaussian random field), the remaining 255 volumes were generated according to the AR1 model X_i=ρ X_i-1+E_i, i= 1, …, 255, where | ρ|< 1 represents a set of random numbers sampled in the 128×128×64 grid and E_i is a set of independent and identically distributed Gaussian random fields. Finally, an additional temporal smoothing was added to these three types of datasets. This was achieved by a temporal convolution with a Gaussian kernel of temporal FWHM equal to 4, 8, and 12 sec, respectively.

Results

PCA of resting-state fMRI data

PCA was applied separately to each single fMRI run (68 runs=17 subjects ×4 runs) and the correlation coefficient between each eigenvariate and the GAS was calculated.

Figure 1A presents the grand average (across the 68 data-sets) of the proportion of the total variance explained by each PC. Figure 1B shows a histogram plot corresponding to the maximum (across eigenvariates within a scan) of the correlation coefficients between the GAS and each of the eigenvariates. The mean value of the maximal correlation coefficient between the GAS and a single PC was r=0.97±0.05 (mean computed over 68 runs), with the smallest maximal correlation being r=0.72. Figure 1C presents the PC index where the maximum value was reached. The first eigenvariate showed maximal correlation with the GAS in 66 out of 68 runs. The second component showed the maximal correlation in only 2 runs (Fig. 1C; below we use the term “atypical” when referring to these 2 runs). The findings presented in Figure 1B and 1C corroborate our hypothesis that PCA is capable of decomposing the data into noncorrelated global effect and remaining network-specific resting-state fluctuations.

FIG. 1.

(A) Grand average (across 68 samples) of the proportion of total variance explained by each PC. (B) Histogram plot of the maximal correlation coefficient between the GAS and each PC (one maximal coefficient is presented for each of the 68 runs). (C) Histogram of the index of the PC that correlated best with the GAS. (D) The grand average of proportion of total variance accounted for by the GAS and the PC-based estimator of the global effect. (E) Grand average (across 68 samples) of the power spectra corresponding to the GAS and the PCA-based global effect estimator. GAS, global average signal; PCA, principal component analysis; PC, principal component.

Next, we computed the proportion of the total variance explained by the GAS and the PC-based global effect estimator (Fig. 1D; mean over all subjects and scans). By definition of PCA, the 1st PC accounts for more variance than any other linear combination of the PCs. Indeed, as expected, our proposed PC-based global effect estimator accounted on average for a larger proportion of the total variance (21.56%±7.32%) than the GAS did (21.06%±7.46%), although this difference was not statistically significant.

The grand average (over the 68 datasets) of the power spectrum corresponding to the GAS and the PC-based global effect estimator indicates that both signals possess similar spectral behavior, reaching the maximum power spectrum value at the frequency of 0.0146 Hz (Fig. 1E).

As seen in Figure 1, in 66 out of the 68 runs the highest correlation coefficient between the GAS and a single PC was obtained between the GAS and the 1st PC. Only 2 out of 68 runs produced an atypical behavior in the sense that the maximal correlation coefficient between the GAS and a single PC was attained with the 2nd PC rather than with the 1st PC. To investigate the differences between those 2 atypical runs and the remaining 66 typical runs, we split the data-set into the two sub-sets. One sub-set included the typical 66 runs in which the first PC showed the highest (compared to all other PCs) correlation with the GAS. The other sub-set included only the above-mentioned 2 atypical cases. A more detailed view of the proportion of the total variance explained by the GAS and each of the two first PCs for these subsets is shown in Table 1.

Table 1.

Proportion of Total Variance Explained by the GAS and the PC-Based Estimator

			% Total variance explained
Subset	Correlation (GAS, PC estimator)	Correlation (GAS, RS fluctuations)	GAS	PC estimator	1st PC	2nd PC
68	0.97±0.05	0.18±0.12	21.06±7.46	21.56±7.32	21.62±7.21	8.15±1.69
Typical 1–66	0.98±0.04	0.17±0.10	21.46±7.19	21.98±7.01	21.98±7.01	8.16±1.67
Atypical 1	0.96	0.29	8.80	10.10	11.69	10.10
Atypical 2	0.76	0.65	5.85	5.60	8.24	5.60

Each row corresponds to a subset of the entire data-set, as indicated by the number of included runs on the left most column. The subsets with 68 and 66 runs correspond to the entire data-set (of 68 runs) and to the subset of runs in which the GAS correlated best with the 1st PC, respectively. The runs Atypical 1 and Atypical 2 correspond to the 2 atypical cases where the GAS was correlated best with the 2nd PC. The 2 next columns show the mean and standard deviation of the correlation coefficient between the GAS and the PC-based global effect estimator and between the GAS and the resting-state fluctuations remaining after removing the PC-based global effect estimator, respectively. The columns to the right present the mean and standard deviation of the proportion of the total variance in the fMRI time-series accounted for by the GAS, by the PC-based global effect estimator, and by each of the two first PCs.

Bold font under '1st PC' or '2nd PC' marks the PC used as the PC-based global effect estimator.

GAS, global average signal; PC, principal component; RS, resting-state; fMRI, functional magnetic resonance imaging.

Table 1 shows, in addition, the correlation coefficients between the GAS and the PC-based estimator, as well as between the GAS and the remaining less global, system-specific fluctuations (after regression of the PC-based estimator). As seen in the table, the run labeled “atypical 1” was indeed “atypical” because its GAS was best correlated to the 2nd PC. However, the correlation coefficient between the GAS and the 2nd PC was high (r=0.96), indicating that in that case too, there was one eigenvariate that closely resembled the GAS. In addition, for run “atypical 1” the PC-based estimator (2nd PC) explained more proportion of the total variance (10.10%) than the GAS (8.80%). In contrast, the second atypical run showed a relatively low correlation coefficient (r=0.76) between the GAS and the PC-based estimator (2nd PC) compared to the 66 typical cases (r=0.98±0.04).

A common feature of the 2 atypical runs is that the GAS accounted for a relatively low proportion of the total variance (8.80% and 5.85% for atypical 1 and 2, respectively) compared to the set of remaining 66 typical cases (21.46%±7.19%). Interestingly, even the 1st PC of those atypical cases also explained a relatively low proportion of the total variance. These findings indicate that for those 2 atypical cases, the global effects are not as global as in the 66 typical cases, in the sense that they might not be spatially extended over the brain.

Figure 2 presents the results of PCA for five different runs, the first three of were typical; that is, their GAS showed the maximal correlation with the 1st PC. The first run (Fig. 2A) corresponds to that typical case where the correlation coefficient between the GAS and the 1st PC reached the maximal value (r=0.99) over the sample of 68 runs. Here, the first PC and the GAS explained 34.20% and 34.05% of the total variance, respectively. The time courses corresponding to the GAS (red) and the first eigenvariate (blue) were found to be almost identical. In contrast, the mean time course computed over all remaining components (green) after regressing-out the 1st PC, namely, the less global, system-specific resting-state fluctuations seemed to represent a different process (correlation with the GAS, r=0.07). Figure 2B and C shows two other typical cases corresponding to correlation coefficients between the GAS and the first PC of r=0.95 and r=0.72 (median and lowest over the sample), respectively. Figure 2A and B is similar in the sense that the GAS and the first PC show very similar changes over time. In contrast, in Figure 2C (r=0.72), the GAS and the first PC are less similar. In fact, the correlation coefficient between the GAS and the remaining resting-state fluctuations produced a value comparable to that between the GAS and the first PC (r=0.7). These findings show a trend in which a decreasing correlation between the GAS and the PC-based estimator results in an increase in the correlation between the GAS and the remaining (after regression of the PC-based estimator) resting-state fluctuations. Panels 2D and 2E show the results of the time courses obtained from the 2 atypical runs. The behavior of the signals in Figure 2D is similar to the one observed in Figure 2A and B (r=0.96), while the atypical case presented in Figure 2E (r=0.76) is more similar to the case presented Figure 2C.

FIG. 2.

PCA analysis performed on five different single fMRI runs. The three time-courses in each panel correspond to the GAS, the PC-based global effect estimator, and the average of the remaining specific resting-state fluctuations (after regressing out the PC-based estimator). (A) A typical run where the GAS correlated best with the 1st PC, with a highest correlation value obtained over the sample (r=0.99). (B) A typical run where the GAS correlated best with the 1st PC with the median correlation value over the sample (r=0.95). (C) A typical run in which the GAS correlated best with the 1st PC; however, the correlation was the minimal over the sample (r=0.72). (D) and (E) show the time-courses associated with the two atypical runs, where the GAS correlated best with the 2nd PC. fMRI, functional magnetic resonance imaging.

Figure 3A presents the t SPMs of the spatial contributions to the two global effect estimators obtained from a hierarchical analysis of all 68 runs from all subjects within our sample. As previously described by Fox and associates (2009), the regions that contribute to the GAS are widespread, extending over large parts of the gray matter. Our proposed PCA-based global effect estimator shows a spatial pattern similar to that presented by the GAS. Figure 3B presents the statistically significant contributions (multiple-comparisons adjusted threshold at α=0.05) to the two global effect estimators superimposed on T1 anatomical images. For both global effect estimators, the ventricles showed up as the regions with the smallest contributions to the global effect.

FIG. 3.

(A) Group-averaged t-SPMs presenting the spatial contributions to the GAS and the PC-based global effect estimator (PC showing the highest correlation with the GAS). The spatial contributions to both global effects estimators extend over large parts of the brain covering most of the gray matter. The smallest contributions can be found in the ventricles. The bottom-most panel presents the differences in contributions to the PC-based estimator and to the GAS. (B) Statistically significant contributions to the global effect estimators following the application of a multiple comparisons corrected threshold of t=± 4.65. The two spatial patterns are very similar, with no regions showing statistically significant differences. SPM, Statistical Parametric Map.

The bottom panels of Figure 3 present SPM maps of differences between the spatial distributions of contributions to the two global effect estimators. Low t values were observed over the larger part of the brain. In fact, setting a threshold of t=± 4.65 (with multiple comparisons correction at α=0.05) produced no regions with statistically significant differences. We therefore concluded that our proposed PC-based global effect estimator resembles the GAS not only temporally but also spatially.

To investigate the mechanisms causing the results in the 2 atypical runs to be different than in the 66 typical runs, we analyzed the spatial contributions to their global effect estimators. Recall that in these atypical cases the GAS correlated best with the 2ed PC, and explained a relatively low proportion of the total variance as compared to the rest of the sample. Thus, the contribution to both global effects estimators is expected to be less extended compared to those presented in Figure 3 for the entire sample. Indeed, since PCA is based on detecting signals of maximum variance, poor spatially extended effects are unlikely to correspond to the 1st PC. Figure 4 shows a t SPM of the spatial contributions corresponding to the GAS, the 1st PC, and the 2nd PC for the case “atypical 1.” As expected, the spatial contributions to the GAS (top part of the Fig. 4) and the PC-based global effect estimator (the 2nd PC here, at the bottom part of the Fig. 4) cover a smaller part of the brain compared to those presented for the entire sample in Figure 3.

FIG. 4.

(A) t-SPMs presenting the spatial contributions to the GAS, the 1st, and the 2nd PC corresponding to the atypical case 1. For this case, the 1st PC explains a relatively low proportion of the total variance as compared to the rest of the sample. (B) Statistically significant contributions to the global effect. The contribution to the GAS and the 2nd PC (PC-based estimator) is less spatially extended compared to those of the typical cases. The 1st PC captures the interaction between the Default Mode Network and the Task Positive Network.

Moreover, in atypical cases the global effect may show a lower variance than the variance associated with a specific network that extends over widespread cortical regions. That network may then explain more of the variance than the GAS and the PC-based global effect estimator. Indeed, Figure 4B shows that after accounting for multiple comparisons with a threshold of t=± 4.65, the 1st PC captured the interaction between the Default Mode Network (Task Negative Network) and the Task Positive Network. A similar pattern of spatial contributions to the 1st PC was obtained also for atypical run 2. We concluded that for both atypical runs, the variance accounted for by the GAS was low, and the 1st PC captured the interaction between the Default Mode Network and the Task Positive Network.

Numerical simulations

The aim of this section is to test whether our finding of high-correlation between the GAS and the 1st PC is a mathematical necessity associated with PCA, or alternatively, a feature intrinsic to the spatio-temporal nature of resting-state BOLD fluctuations. We therefore test the hypothesis that when pursuing PCA, there must be one eigenvariate highly correlated with the GAS. We also tested whether regressing out the eigenvariate that correlates with the GAS with the highest correlation value compared to all other eigenvariates mandates spurious negative seed-based correlations, as regression of the GAS does.

For each of the four types of simulated datasets (temporally uncorrelated and temporally smoothed AR1 models), we generated a set of 100 independent samples, each consisting of a time series of 256 volumes. PCA was carried out on each individual sample. Then, the GAS (spatial average over the 128×128×64 grid) was compared with each of the eigenvariates produced by the PCA. Each row in Figure 5 corresponds to a different type of simulated dataset. For each type, the first column shows the grand average (across the 100 samples) of the proportion of the total variance explained by each PC. Note that, similarly to the measured fMRI dataset, most of the total variance is accounted for by the first 50 components, although this number seems to decrease with increasing temporal smoothing. The second column presents a histogram plot corresponding to the maximum (across eigenvariates) of the correlation coefficients between the GAS and each of the eigenvariates. The third column presents a bar plot for the PC index in which the maximum value of these correlation coefficients was reached. Note that the maximal correlation values between the simulated GAS and single PCs were significantly lower than those obtained for the fMRI data (compare the panels in the second column of Fig. 5 to Fig. 1B). In addition, when comparing the indices within the third column of Figure 5 to those in Figure 1C, the maximum of these correlation coefficients were not necessarily obtained in the first two PCs. With these results, we concluded that when pursuing PCA, there must not necessarily be one eigenvariate highly correlated with the GAS. Therefore, the high correlation observed in resting-state fMRI data between the first PC and the GAS was not an effect mandated by PCA; rather, it is a feature of fMRI-measured resting-state activity.

FIG. 5.

PCA applied to four different types of simulated datasets. The first row corresponds to temporally uncorrelated volumes. The next three rows present results of datasets simulated with temporal autocorrelation following AR1 model and additional temporal smoothing at three different scales. For each case, the first column shows the grand average (over 100 datasets) of the proportion of the total variance explained by each PC. Note that the proportion of the total variance explained by the first PC increases with the temporal smoothing of the data. The second column shows the maximal correlation coefficient attained between the GAS and each PC. The third column presents the PC indices at which such maximal values were reached. Note that the values of these maximal correlations are not as high as those obtained with actual resting-state fMRI data presented in Figure 1. The fourth column shows histograms of the spatial sum of the coefficients β_f in the PC-based estimator regression. In contrast to the removal of the GAS, the sum of these coefficients is evidently different from zero.

Finally, it has been demonstrated (Fox et al., 2009; Murphy et al., 2009) that the appearance of negative correlations following the removal of the GAS is associated with the outcome that the spatial sum of the coefficients β_f in model (2) is equal to zero. To verify that regressing-out the PC-based estimator does not mandate spurious negative correlations, we selected a seed time course corresponding to the voxel with coordinates [64 64 32]. We then computed the correlation coefficients between the simulated time-course of that point and all other points in the grid after regressing-out the PC-based estimator according to model (2). The fourth column in Figure 5 shows histograms (computed over the 100 simulated runs) of the spatial sum of the coefficients β_f , integrated over the grid. In contrast to the removal of the GAS, the sum of the β_f coefficients here is evidently different from zero, which is an indicator that the removal of the PC-based estimator does not necessarily mandate negative correlations.

Global effects and anti-correlations

Seed correlation maps were generated after separately using the two strategies for removal of global effect: regressing out the GAS and the proposed PC-based global effect estimator. Figures 6 and 7 show axial-oblique slices of these correlation maps with a posterior cingulate cortex (PCC) seed and an MT/V5 seed, respectively. These figures present the correlations obtained using the data from all subjects within our sample, by generating a hierarchical random effect model. The panels to the right present the regions showing statistically significant positive (in red) and negative (in blue) correlations with the seed time-course following multiple comparisons corrections (with a threshold at t=5.17). In both Figures 6 and 7, we observed statistically significant anti-correlations between areas of the Default Mode Network (DMN) and the Task Positive Network (TPN).

FIG. 6.

(A) Group-averaged correlation maps (using a t-SPM scale) of voxels showing correlation with the PCC seed. (B) Statistically significant anti-correlations were observed between the Task-Positive Network and the Default Mode Network by applying a multiple comparisons corrected threshold of t=5.17 following the regression of either the GAS or the PC-based global effect estimator. The anti-correlations obtained following the regression of the PC-based global effect estimator were less extended than those obtained by regression the GAS. PCC, posterior cingulate cortex.

FIG. 7.

(A) Group-averaged correlation maps (using a t-SPM scale) of voxels showing correlation with the middle temporal (MT) seed. (B) Statistically significant anti-correlations were observed between the Task-Positive Network and the Default Mode Network by applying a multiple comparisons corrected threshold of t=5.17 following the regression of either the GAS or the PC-based global effect estimator. As in Figure 6, the anti-correlations obtained following the regression of the PC-based global effect estimator were less extended than those obtained by regression the GAS.

For the two seed ROIs used, anti-correlations between the default mode and the task-positive networks appeared not only following the regression of the GAS but also following the regression of the PCA-based global effect estimator. These findings confirm previous reports that the DMN and the TPN show anti-correlated fluctuations in the resting state. Note that the negative correlations obtained following the removal of the PC-based estimator were less spatially extended then those appearing following the regression of the GAS. Given that no statistically significant differences between the spatial contributions to the two global effect estimators were observed (Fig. 3), one could expect that regressing out the PC-based estimator might cause spurious negative correlations, just as regressing out the GAS does. Similar expectations could be raised because there were virtually no statistically significant differences between the correlation maps obtained for the MT seed after regressing out the GAS and the PC-based estimator (Fig. 7). However, for that MT seed case, the spatial sum of the coefficients β_f in model (2) after regressing-out the PC-based estimator was 1205.8±4615.5 (mean±SD computed over 68 runs; percentile 5%, −6836.2; percentile 95%, 8055.6). The corresponding spatial sum for the PCC seed case was 2670.0±5149.7 (percentile 5%, −6384.9; percentile 95%, 10648.5). In both cases, the spatial sum of the coefficients β_f was different than 0, indicating (but not confirming, please see details in the Discussion section) that the PC-based estimator does not impose spurious negative correlations. In contrast, the sum of the coefficients β_f after regressing out the GAS was −2.59e-007±3.26e-006 (percentile 5%, −5.38e-006; percentile 95%, 6.76e-006) for the MT seed. The corresponding sum for the PCC seed was −3.71e-007±3.24e-006 (percentile 5%, −4.65e-006; percentile 95%, 7.05e-006). Note that after regressing-out the GAS we obtained an approximate rather than exact zero value for the sum of the coefficients β_f . The reason is that our seed points were not defined as single voxels but as the average of voxels within 6 mm spherical regions around predetermined locations. In both cases, these approximate zero values confirmed that regressing out the GAS might introduce spurious negative correlations.

These findings corroborate our hypothesis that regression of the GAS is a source of spurious anti-correlations between the Default Mode Network and the Task Positive Network beyond those allowed by the RSLAM, while the removal of the PC-based estimator does not mandate spurious negative correlations.

Discussion

Summary of the results

The results from our study show that PCA decomposition of resting-state fMRI data produces a single temporal component closely resembling the temporal pattern of the GAS. Therefore, this eigenvariate can be considered a good estimator of global effects. Due to the basic properties of PCA, by definition this eigenvariate is uncorrelated with the remaining PCs that represent the more network-specific resting-state fluctuations. In addition, we demonstrated that regressing out the eigenvariate that closely resembles the GAS produces seed-based negative correlations between the default mode network and the task-positive network. Since our PC-based estimator does not introduce spurious correlations at the estimation level, these negative correlations do exist, rather than reflecting an analysis-induced artifact.

An eigenvariate closely resembling the GAS is not a mathematical necessity

By means of numerical simulations we have shown that having a single PC in general, and the first PC in particular, closely resembling the GAS is not a mathematical necessity of PCA (Fig. 5, columns 1–3). Our simulations are not meant to mimic BOLD resting-state fluctuations. Instead, they are meant to act as a counterexample, refuting the hypothesis “there is always an eigenvariate with a close, standing-out resemblance to the GAS compared to all other PCs, as mandated by the algebraic manipulations involved in the PCA, independent of the type of signal.” In other words, the simulations are meant to refute a claim related to PCA, not necessarily related to resting-state BOLD-fMRI. We therefore applied a simple scenario, not necessarily mimicking BOLD-fMRI resting-state fluctuations in a precise manner, but sufficient for refuting the aforementioned hypothesis.

These simulations demonstrate that our results do not merely reflect an effect derived from the algebraic manipulations involved in PCA; rather, they reflect an intrinsic property of the complex spatio-temporal nature of the fMRI-measured resting-state data.

Global and network-specific fluctuations are uncorrelated, supporting an RSLAM

Removal of the GAS in resting-state studies has been based on an implicitly hypothesized additive model. This model assumes that the fMRI-measured resting-state fluctuations are given by the super-position of a global effect signal and fluctuations that reflect specific interactions of major neuronal systems. A first step to support such a hypothesis was pursued by Fox and associates (2009), who showed that the GAS is not simply the average of fluctuations with origin in major brain networks. However, the degree of correlation between the GAS and the specific activity of resting-state networks has not been elucidated, and the hypothesized additive model has not been properly evaluated (Murphy et al., 2009).

In the current article, we propose a PC-based estimator of the global effect in fMRI resting-state data, better suited for the role than the currently used GAS estimator. Specifically, our results show that PCA decomposes fMRI-measured resting-state activity into a component that estimates the global effect and the remaining components, whose sum approximates the lumped system-specific resting-state fluctuations. The fit of the PC-based estimation to the data (in 66/68 cases the 1st PC correlated best with GAS, with high correlation values, mean=0.97±0.05) makes a strong case that PCA is an appropriate method for identifying the global effect and separating it from the remaining network-specific signals. The two components are not correlated because they consist of different, orthogonal PCs. Moreover, as discussed above, our simulations show that these features of the global effect are inherent to the spatiotemporal dynamics of resting-state fMRI data, rather than being mandated by PCA.

It is worth noting that there is no ground truth on which to rely when considering what portions of the signal are global and what portions are system specific. In fact, the precise determination of the global component and the remaining system-specific components is subject to definition. However, the similarity of the PC-based estimator to the GAS and the fact that it is uncorrelated with the remaining network-specific signals make it a natural, data-driven candidate for this role.

The fit of this PCA-based model to the data, showing that a single PC approximately captures the global effect, and the fact that the global effect is uncorrelated with the remaining network-specific fluctuations support the concept of the RSLAM we formally defined. Overall, our findings support an RSLAM model in which a global effect and system-specific fluctuations are super-positioned and uncorrelated.

PC-based global effect estimator, RSLAM, and the origin of the global effect

We emphasize that the concept of the RSLAM and regressing out the PC-based global effect estimator are proposed independently of whether neural activity-related resting-state signals contribute to the global signal. Neither regression of the GAS nor PCA can shed light on what comprises the global signal. We propose a model of uncorrelated global effect and system-specific components that fits resting-state fMRI data very well. However, we posit that this fit to the data is independent of whether the global signal reflects neurophysiological activity or systemic physiological effects or both. The method we propose for removing the global effect signal can be used for facilitating the observation of network-specific functional connectivity in both cases, when the global signal does or does not reflect neurophysiological activity.

PC-based resting-state global effect estimator does not assign each network to a PC

We would like to emphasize that we do not claim to assign individual networks to single PCs. Our only use of PCA is to decompose the signal to a global effect and all remaining signals lumped together. Indeed, Beckmann and associates (2005) showed that the mathematical constraint of orthogonality within the set of spatial components (eigen images) does not necessarily imply that large areas of activation that overlap significantly with other spatial components can no longer be extracted. This is exactly what we observe in our results: a PC highly similar to the GAS could be extracted from the data even by imposing (due to PCA) the restrictive condition of orthogonality between spatial components.

PC-based resting-state global effect estimator: comparison to previous task-based studies

Our PCA-based global effect estimator is related to the model and the removal strategy proposed by Macey and associates (2004) for response-based and hypercapnia studies. Although Macey's method is not specifically designed for fMRI resting state and does not rely on PCA, it is based on the assumption that a global effect is replicated with the same temporal pattern throughout the brain, although not necessarily at the same magnitude. Indeed, our analysis shows that the time course of the eigenvariate that correlates best with the GAS is replicated across the brain at different magnitudes, which are determined by the values of the corresponding eigenimage. The strategy taken here is related, in addition, to the global signal removal method presented by Andersson and associates (2001) for PET activation studies. These authors showed that rotation matrices applied to the eigenvariates obtained by PCA produced conceivable candidates for representing the global effect (Andersson et al., 2001). The method was based on maximizing the spatial contributions to the global effect while discounting contributions from any variable with a highly localized response (i.e., restricting the stimulus-dependent activations to be very local and limited to a small number of distinct foci). However, these assumptions cannot be guaranteed in fMRI resting state since networks-specific fluctuations extend over wide-spread functional systems. Here we showed that a single PC captures the global effect in resting-state fMRI.

On the possibility of an independent component analysis-based global effect estimator

An alternative decomposition that may fit the global effect could use independent component analysis (ICA). The advantage of ICA is that, unlike PCA, it does not enforce orthogonality between the estimated global effect and the sum of the network-specific fluctuations. PCA is more restrictive than ICA in the sense that it enforces orthogonality of the components in both the spatial and temporal domains. To adapt our methodology to ICA, we would have to compare the GAS with temporal ICs and select that single IC that is (i) highly correlated with the GAS and (ii) accounts for a large proportion of the total variance, while (iii) keeping independence among all ICs. The latter condition is necessary to avoid introduction of artificial correlations, which makes temporal ICA type method a better choice compared to spatial ICA. Another reason why temporal ICA could be a good option is that the effectiveness of spatial ICA depends more on the sparseness of the extracted components than on their independence (Daubechies et al., 2009). In other words, spatial ICA has been proven to be effective at detecting sparse signals, and clearly a global-effect component is, by definition, not spatially sparse. It follows that if an ICA is to be done, a temporal ICA technique should be selected, because spatial ICA would be ineffective. Note, however, that any candidate temporal IC fulfilling (i), (ii), and (iii) would express less proportion of the total variance than the 1st PC does, since any single IC can be expressed as a linear combination of all PCs, which by definition of PCA explains less variance than the 1st PC. Therefore, we would need to establish a tradeoff between conditions (i) and (ii). That is, we could obtain a single IC more correlated to the GAS than our PC-based estimator, but this would come at the expense of decreasing the proportion of the total variance explained (as compared to the PC-based estimator). It would be of interest to check whether the estimation of the global effect could be improved while establishing a tradeoff between conditions (i) and (ii) and the use of either temporal ICA or Anderson's rotation of PCA matrix (Andersson et al., 2001).

Sources of modified seed-based correlations introduced by regressing out the global effect

The goal of seed-based functional connectivity analysis is to detect regions in the brain that are functionally connected with the seed region. Regressing out the GAS in resting-state studies has enhanced the network-specific seed-based correlations (Fox et al., 2009). The GAS has been the common choice for a global effect estimator in such a regression strategy. However, as pointed out by Murphy and associates (2009), a major drawback is that spurious negative correlations might be artificially introduced by the regression of the GAS. By definition, the GAS is correlated with the remaining less-global, network-specific components, since it mixes signals from all voxels and networks. This is a critical source of artificial negative correlations in seed-based functional connectivity analysis (Murphy et al., 2009). Here we provide a more detailed view on the origin of the spurious negative correlations detected by Murphy and associates (2009). From the statistical point of view, there are two causes to the modified correlations values following the regression of the GAS: modeling and estimation. With regard to the modeling level, regression of the GAS implicitly relies on an assumed RSLAM. This model assumes that the fMRI-measured resting-state fluctuations are given by a global effect signal added on top of the system-specific fluctuations. Regression of the global effect decreases the overall seed-based correlation values to a level justified by the RSLAM.

With regard to the estimation level, a second source of spurious correlations arises when estimating the global effect as the GAS. This causes an additional decrease in correlation values at the estimation level, beyond what is justified by the RSLAM. In contrast, our PC-based global effect estimator does not introduce artificial correlations beyond the decrease in correlation values allowed by the RSLAM. Our PC-based estimator meets the requirements imposed by the RSLAM, since this estimator is by definition orthogonal to the remaining system-specific fluctuations, avoiding the spurious correlations introduced by the regression of the GAS. In other words, the PC-based estimator decreases correlation values only up to the level justified by the RSLAM, avoiding any additional modification of correlation values due to estimation. Therefore, our proposed PC-based global effect estimator is better suited as a global effect estimator than the GAS is.

Regressing out the GAS from the original fMRI time-series sets the sum of beta coefficients (β_f ) from equation (2) to zero. This condition mandates the occurrence of spurious negative correlations. Note, however, that a sum of beta coefficients different from zero after regressing out a global effect estimator does not necessarily guarantee the avoidance of spurious negative correlations. Avoiding spurious negative correlations introduced at the estimation level is guaranteed by regressing out a global effect estimator that is uncorrelated with the remaining fluctuations. Therefore, the artifact-free regression of our PC-based estimator relies on the fact that it is uncorrelated with the remaining PCs.

Although we demonstrate similar spatial distributions of contributions to the GAS and the PC-based estimator, we strongly recommend the use of the PC-based global effect estimator rather than the GAS. The use of the PC-based estimator avoids the theoretical limitations imposed by using the GAS and potentially misleading interpretations of functional connectivity results.

Default-mode and task-positive networks are anti-correlated

We have shown that using the PC-based estimator of the global effect confirms the existence of negative correlations between the network-specific resting-state activities of regions in the Default Mode and the Task Positive networks. Importantly, based on the RSLAM, regression of our proposed PC-based global effect estimator eliminates a source of artificial negative correlations related to the estimation level. Our results support the concept that these negative correlations are intrinsic to the brain, rather than merely being an artifact introduced by the regression of the global effect. In fact, resting-state anti-correlated networks have been already detected without the regression of the GAS (Chai et al., 2012; Chang and Glover, 2009; Fox et al., 2009). In addition, these two networks show anti-correlated signals during task-based paradigms. Thus, this finding adds yet another feature of similarity (anti-correlations in task-evoked paradigms as well as in resting-state) to the already-reported similarity between task-evoked networks and resting-state networks (Smith et al., 2009).

Conclusions

We have used PCA to define a novel estimator of the global effect in resting-state fMRI studies. Our findings demonstrate that the global effect can be captured by a single PC, which, by definition, is uncorrelated with the remaining resting-state components. Our results therefore support the hypothesis of the RSLAM in which global effect and system-specific fluctuations are super-positioned and uncorrelated. Our proposed estimator of the global effect can be regressed out from resting-state fMRI time-series without introducing spurious negative correlations beyond the global correlation decrease allowed by the assumed linear additive model. Using our PCA-based method, we have shown that the previously reported anti-correlation between fluctuations in the task-positive and default mode networks is an intrinsic property of the brain, rather than an artifact induced by the regression of the GAS.

Footnotes

Acknowledgments

Supported by CIHR Grant MOP-102599 and Human Frontier Science Grant RGY0080/2008 awarded to A.S., and by Industry Canada/MNI Center of excellence in commercialization and research postdoctoral fellowship and grant awarded to F.C. and A.S., respectively. We thank Debra Dawson and Laura Betcherman for their comments on the article and English editing.

Authors' Comment

An implementation of the PC-based estimator of the global effect and its regression is available in NIAK (www.nitrc.org/projects/niak).

Author Disclosure Statement

No competing financial interests exist.

*

http://code.google.com/p/niak/

†

www.math.mcgill.ca/keith/fmristat/

References

Aguirre

, Zarahn

, D'Esposito

. 1998. The inferential impact of global signal covariates in functional neuroimaging analyses. Neuroimage, 8:302–306.

Andersson

, Ashburner

, Friston

. 2001. A global estimator unbiased by local changes. Neuroimage, 13:1193–1206.

Baumgartner

, Ryner

, Richter

, Summers

, Jarmasz

, Somorjai

. 2000. Comparison of two exploratory data analysis methods for fMRI: fuzzy clustering vs. principal component analysis. Magn Reson Imaging, 18:89–94.

Beall

, Lowe

. 2007. Isolating physiologic noise sources with independently determined spatial measures. Neuroimage, 37:1286–1300.

Beckmann

, DeLuca

, Devlin

, Smith

. 2005. Investigations into resting-state connectivity using independent component analysis. Philos Trans R Soc B, 360:1001–1013.

Bianciardi

, Fukunaga

, van Gelderen

, Horovitz

, de Zwart

, Shmueli

, Duyn

. 2009. Sources of functional magnetic resonance imaging signal fluctuations in the human brain at rest: a 7 T study. Magn Reson Imaging, 27:1019–1029.

Birn

, Diamond

, Smith

, Bandettini

. 2006. Separating respiratory-variation-related fluctuations from neuronal-activity-related fluctuations in fMRI. Neuroimage, 31:1536–1548.

Biswal

, Yetkin

, Haughton

, Hyde

. 1995. Functional connectivity in the motor cortex of resting human brain using echo-planar MRI. Magn Reson Med, 34:537–541.

Carbonell

, Bellec

, Shmuel

. Artifact-Free Removal of Global Signal in Resting-State Reveals Anti-Correlated Networks. 16th Annual Meeting of the Organization of the Human Brain Mapping, Barcelona, Spain, 2010program number 975.

10.

Chai

, Castanon

, Ongur

, Whitfield-Gabrieli

. 2012. Anticorrelations in resting state networks without global signal regression. Neuroimage, 59:1420–1428.

11.

Chang

, Glover

. 2009. Effects of model-based physiological noise correction on default mode network anti-correlations and correlations. Neuroimage, 47:1448–1459.

12.

Collins

, Neelin

, Peters

, Evans

. 1994. Automatic 3D intersubject registration of MR volumetric data in standardized Talairach space. J Comput Assist Tomogr, 18:192–205.

13.

Daubechies

, Roussos

, Takerkart

, Benharrosh

, Golden

, D'Ardenne

, Richter

, Cohen

, Haxby

. 2009. Independent component analysis for brain fMRI does not select for independence. Proc Natl Acad Sci, 106:10415–10422.

14.

Desjardins

, Kiehl

, Liddle

. 2001. Removal of confounding effects of global signal in functional MRI analyses. Neuroimage, 13:751–758.

15.

Fox

, Raichle

. 2007. Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging. Nat Rev Neurosci, 8:700–711.

16.

Fox

, Snyder

, Vincent

, Corbetta

, Van Essen

, Raichle

. 2005. The human brain is intrinsically organized into dynamic, anticorrelated functional networks. Proc Natl Acad Sci, 102:9673–9678.

17.

Fox

, Zhang

, Snyder

, Raichle

. 2009. The global signal and observed anticorrelated resting state brain networks. J Neurophysiol, 101:3270–3283.

18.

Frank

. 2000. Impact of a confounding variable on a regression coefficient. Sociol Methods Res, 29:147–194.

19.

Fransson

. 2005. Spontaneous low-frequency BOLD signal fluctuations: an fMRI investigation of the resting-state default mode of brain function hypothesis. Hum Brain Mapp, 26:15–29.

20.

Friston

, Frith

, Liddle

, Dolan

, Lammertsma

, Frackowiak

. 1990. The relationship between global and local changes in PET scans. J Cereb Blood Flow Metab, 10:458–466.

21.

Friston

, Frith

, Liddle

, Frackowiak

. 1993. Functional connectivity: the principal-component analysis of large (PET) data sets. J Cereb Blood Flow Metab, 13:5–14.

22.

Gavrilescu

, Shaw

, Stuart

, Eckersley

, Svalbe

, Egan

. 2002. Simulation of the effects of global normalization procedures in functional MRI. Neuroimage, 17:532–542.

23.

Greicius

, Krasnow

, Reiss

, Menon

. 2003. Functional connectivity in the resting brain: a network analysis of the default mode hypothesis. Proc Natl Acad Sci, 100:253–258.

24.

Lowe

, Mock

, Sorenson

. 1998. Functional connectivity in single and multislice echoplanar imaging using resting-state fluctuations. Neuroimage, 7:119–132.

25.

Lund

, Madsen

, Sidaros

, Luo

, Nichols

. 2006. Non-white noise in fMRI: does modelling have an impact? Neuroimage, 29:54–66.

26.

Macey

, Macey

, Kumar

, Harper

. 2004. A method for removal of global effects from fMRI time series. Neuroimage, 22:360–366.

27.

Murphy

, Birn

, Handwerker

, Jones

, Bandettini

. 2009. The impact of global signal regression on resting state correlations: are anti-correlated networks introduced? Neuroimage, 44:893–905.

28.

Otsu

. 1979. A thresholding selection method from gray-level histogram. IEEE Trans Syst Man Cybern, 9:62–66.

29.

Perlbarg

, Bellec

, Anton

, Pelegrini-Issac

, Doyon

, Benali

. 2007. CORSICA: correction of structured noise in fMRI by automatic identification of ICA components. Magn Reson Imaging, 25:35–46.

30.

Schölvinck

, Maier

, Ye

, Duyn

, Leopold

. 2010. Neural basis of global resting-state fMRI activity. Proc Natl Acad Sci, 107:10238.

31.

Shmuel

, Leopold

. 2008. Neuronal correlates of spontaneous fluctuations in fMRI signals in monkey visual cortex: implications for functional connectivity at rest. Hum Brain Mapp, 29:751–761.

32.

Smith

, Fox

, Miller

, Glahn

, Fox

, Mackay

, Filippini

, Watkins

, Toro

, Laird

, Beckmann

. 2009. Correspondence of the brain's functional architecture during activation and rest. Proc Natl Acad Sci U S A, 106:13040–13045.

33.

Weissenbacher

, Kasess

, Gerstl

, Lanzenberger

, Moser

, Windischberger

. 2009. Correlations and anticorrelations in resting-state functional connectivity MRI: a quantitative comparison of preprocessing strategies. Neuroimage, 47:1408–1416.

34.

Wise

, Ide

, Poulin

, Tracey

. 2004. Resting fluctuations in arterial carbon dioxide induce significant low frequency variations in BOLD signal. Neuroimage, 21:1652–1664.

35.

Worsley

, Cao

, Paus

, Petrides

, Evans

. 1998. Applications of random field theory to functional connectivity. Hum Brain Mapp, 6:364–367.

36.

Worsley

, Liao

, Aston

, Petre

, Duncan

, Morales

, Evans

. 2002. A general statistical analysis for fMRI data. Neuroimage, 15:1–15.

37.

Zarahn

, Aguirre

, D'Esposito

. 1997. Empirical analysis of BOLD fMRI statistics. I. Spatially smoothed data collected under null-hypothesis and experimental conditions. Neuroimage, 5:179–197.

38.

Zijdenbos

, Forghani

, Evans

. 2002. Automatic “pipeline” analysis of 3-D MRI data for clinical trials: application to multiple sclerosis. IEEE Trans Med Imaging, 21:1280–1291.