|Home | About | Journals | Submit | Contact Us | Français|
Low-frequency fluctuations in fMRI signal have been used to map several consistent resting state networks in the brain. Using the posterior cingulate cortex as a seed region, functional connectivity analyses have found not only positive correlations in the default mode network but negative correlations in another resting state network related to attentional processes. The interpretation is that the human brain is intrinsically organized into dynamic, anti-correlated functional networks. Global variations of the BOLD signal are often considered nuisance effects and are commonly removed using a general linear model (GLM) technique. This global signal regression method has been shown to introduce negative activation measures in standard fMRI analyses. The topic of this paper is whether such a correction technique could be the cause of anti-correlated resting state networks in functional connectivity analyses. Here we show that, after global signal regression, correlation values to a seed voxel must sum to a negative value. Simulations also show that small phase differences between regions can lead to spurious negative correlation values. A combination breath holding and visual task demonstrates that the relative phase of global and local signals can affect connectivity measures and that, experimentally, global signal regression leads to bell-shaped correlation value distributions, centred on zero. Finally, analyses of negatively correlated networks in resting state data show that global signal regression is most likely the cause of anti-correlations. These results call into question the interpretation of negatively correlated regions in the brain when using global signal regression as an initial processing step.
Spontaneous fluctuations in blood oxygenation level dependent (BOLD) fMRI signals have recently aroused a large amount of interest in the fMRI literature (Fox and Raichle, 2007). These fluctuations are often correlated between functionally related areas and can occur either on top of task-induced signal modulations (Fox et al., 2006) or in the absence of an explicit task (Biswal et al., 1995). It has been hypothesized that correlated fluctuations reflect synchronized variations in the neuronal activity of discrete brain areas and are representative of functional connections within networks of the brain. Functional connectivity analyses can investigate these coherent signal fluctuations, characterized by their low frequency (~0.1 Hz) (Lowe et al., 1998). By studying fluctuations at rest, researchers have claimed that the brain is intrinsically organized into dynamic, anti-correlated functional networks (Fox et al., 2005; Fransson, 2005; Greicius et al., 2003). The extent to which a commonly used correction method called global signal regression alters functional connectivity maps by introducing anti-correlated time series is the topic of this paper.
Biswal et al. first reported the correlation in fMRI signal fluctuations between the left and right motor cortices when the brain was at “rest” (Biswal et al., 1995). Subsequent studies have identified several consistent resting state networks, including motor, auditory, visual, attention and default mode (Damoiseaux et al., 2006; De Luca et al., 2006; Greicius et al., 2003). The default mode network is of particular interest since it appears to be more active during rest than during task (Raichle et al., 2001). It has been hypothesized that this activation is indicative of internal monitoring, e.g. “day-dreaming” or a “wandering mind” (Buckner et al., 2008; Mason et al., 2007). Early clinical studies have attributed disruptions in the connections between nodes of the default mode network to disorders such as Alzheimer’s disease (Greicius et al., 2004), schizophrenia (Garrity et al., 2007), ADHD (Sonuga-Barke and Castellanos, 2007) and autism (Just et al., 2007; Kennedy et al., 2006).
The common assumption in most fMRI investigations of connectivity is that correlated fluctuations in resting state networks are neuronal in origin. However, other sources of fluctuations exist in fMRI data that are not directly related to local neuronal firing and have a largely physiological source. Cardiac pulsations and respiration-related artefacts can cause significant correlated signal changes in the vicinity of large blood vessels and throughout grey matter, obscuring spontaneous neuronal fluctuations (Lund et al., 2006). By measuring cardiac and respiratory traces during acquisition of fMRI data, the influence of the related physiological fluctuations can largely be removed, using retrospective correction techniques (Glover et al., 2000; Hu et al., 1995). Other techniques have been developed that do not require physiological recordings (Beall and Lowe, 2007; Chuang and Chen, 2001; Perlbarg et al., 2007). However, further physiological sources of noise in the fMRI signal remain uncorrected by these methods. In particular, low frequency BOLD fluctuations have been associated with changes in the level of arterial carbon dioxide (CO2) in the frequency range 0–0.05 Hz (Wise et al., 2004) resulting from changes in the breathed-volume (respiration volume over time — RVT) (Birn et al., 2006) and changes in the rate of cardiac pulsation at ~0.08 Hz (Shmueli et al., 2007). These sources of noise can artificially inflate connectivity measures since they introduce global, spatial coherence across the brain.
Global signal regression, otherwise known as orthogonalization to the global signal, is often performed as a processing step in an attempt to account for several potential sources of physiological noise. In this technique, the global signal, calculated by averaging the time series over all voxels in the brain, is used as a regressor in a general linear model (GLM) to remove the associated variance (Desjardins et al., 2001; Macey et al., 2004). This technique assumes that fMRI experiments are concerned with local changes in neuronal activity and that global signals represent uninteresting sources of noise. However, this assumption is only accurate when the global signal and experimental conditions are orthogonal to each other, that is, uncorrelated. The majority of the methodology literature regarding global signal regression examines its effect on task paradigms. Whether findings from such research are applicable to seed-based connectivity studies is unknown.
The global signal correction technique has been widely used in the PET imaging literature to remove global blood flow fluctuations with some debate as to the appropriate removal method (Andersson et al., 2001; Arndt et al., 1996; Fox etal., 1988; Friston et al., 1990; Ramsay et al., 1993). It has been suggested that defining the global signal as the average over all voxels can introduce a bias since the resulting time course will not be orthogonal to the task-induced activations (Andersson, 1997; Strother et al., 1995). Given that voxels responding to the experimental task are included in the global regressor, this correction technique will significantly influence the results by underestimating true activation levels and by introducing deactivations.
The global signal explains a large proportion of spatial variance in PET imaging (Fox and Mintun, 1989; Friston et al., 1990). Whether removal of this global signal is prudent in fMRI studies of task-related activity is a matter of debate. In spatially smoothed fMRI data, the global signal was shown to be strongly influenced by the performance of a behavioural task but inclusion of this signal as a covariate did not reduce sensitivity (Aguirre et al., 1997). Also, including the global signal covariate in a general linear model reduced spatial coherence and subsequently stabilised false-positive rates (Zarahn et al., 1997). However it has also been demonstrated that if the global signal is strongly correlated with experimental manipulations, considerably different results may be obtained. Aguirre et al. analysed an fMRI dataset in which subjects made button presses to brief intermittent stimuli both with and without covariation for the global signal and found that global signal regression reduced the spatial extent and intensity of positively activated regions whilst introducing spurious negatively activated areas (Aguirre et al., 1998a). Interestingly, these negative regions map out what is now known as the default mode network, a network shown to be anti-correlated with attentional tasks (i.e., deactivated during attention demanding tasks). Although the authors state that these areas “did not have a significant relationship with the task”, global signal regression may have helped bring these true deactivations to light.
Five different methods of global normalization in fMRI were compared by Gavrilescu et al., namely grand mean session scaling, proportional scaling, ANCOVA, a masking method and an orthogonalization method (Gavrilescu et al., 2002). The orthogonalization method, first proposed by Desjardins et al. (2001), was shown to perform better than other methods and is essentially the same as the global signal regression technique. However, it can, nonetheless, decrease the sensitivity of statistical analyses and induce artefactual task deactivations. Using both a high and low arousal emotional task, Junghofer et al. showed that the validity of proportional global signal scaling varied as a function of the emotional arousal of the stimuli (Junghofer et al., 2005). The high arousal condition violated the assumption of orthogonality between the global regressor and the experimental condition leading to a reduction in positively activated areas and “widely distributed artificial deactivations”, none of which made sense in terms of known emotional and default mode areas.
It is clear that global signal regression can cause reductions in sensitivity and can introduce false deactivations in studies of task activation since the assumption of orthogonality can be violated when the experimentally-induced activations contaminate the global signal. By definition, the experimental condition in resting state data is undefined. Although it is known that resting state fluctuations are low-frequency, the exact timing, spatial extent and relative phase between areas are unknown and may vary from session to session. The degree of correlation between the global signal and the resting state fluctuations cannot therefore be determined, and thus global signal regression could lead to spurious results in seed voxel correlation analyses of resting state data.
In a seed voxel analysis, a common method for computing functional connectivity, a time series is chosen in the brain that is hypothesized to represent fluctuations-of-interest. The correlation value between this time series and every other voxel indicates the extent of functional connectivity between the voxels. Using this connectivity analysis, previous studies have claimed that the brain is intrinsically organized into dynamic, anti-correlated networks (Fox et al., 2005; Fransson, 2005; Greicius et al., 2003; Kelly et al., 2008). These studies, however, have all used global signal regression as a pre-processing step. The default mode network is termed the task-negative network since it deactivates during a wide range of tasks. Its anti-correlated network is termed the task-positive network since it represents regions that increase in activity when default mode areas decrease. If the global signal is uncorrelated with the resting state fluctuations as represented by the seed voxel, then the finding of a task-positive network (anti-correlated with the task-negative network) at rest is valid. If, on the other hand, this is not the case, the interpretation that the brain is organized into anti-correlated networks may be brought into question since it is based on results that are contaminated by the global signal regression preprocessing step.
The extent to which global signal regression affects seed voxel functional connectivity analyses was investigated. This paper is organized into four sections (theory, simulations, breathing and visual data, and resting state data), each of which addresses a different aspect of resting state fluctuations and how they are influenced by global signal regression. The theory demonstrates that mathematically, regardless of the characteristics of the resting state fluctuations, global signal regression will always result in a negative mean correlation value during a seed voxel functional connectivity analysis. Three sets of simulations were performed. The first is an empirical demonstration of the theory. The second reveals that if similar resting state fluctuations (represented by a sine wave) are present in all voxels, functional connectivity analyses will determine that approximately half are anti-correlated with a seed after global signal regression. The relative phase of fluctuations in the seed voxel and other regions can bias those regions to become anti-correlated, serving as a possible explanation for the clustering of anti-correlated areas in the brain. The third simulation investigated the behaviour of these anti-correlations as the spatial extent of voxels containing resting state fluctuations increases from zero voxels to the entire brain. The challenge in resting state functional connectivity studies is that the timing and location of spontaneous neuronal activation and global nuisance fluctuations (such as respiration) are unknown. Since it is impossible to systematically alter resting fluctuations, visual task activation was used in the breathing and visual data section to create localisable connectivity maps similar to those generated in resting state correlation studies. A comparatively global fluctuation was introduced into the data using a breath-holding challenge. The relative phase of the two tasks were varied to examine how the influence of global signal regression on connectivity measures is dependent upon the relationship between resting state fluctuations and the global signal, a quantity which, by definition, is unknown. Finally, task-negative and task-positive regions were defined in resting state data and changes in correlation measure with global signal regression were examined in depth.
This mathematical proof demonstrates that after global signal regression, the sum of correlation values with a seed voxel across the entire brain is less than or equal to zero.
Let each Si(t) be a column vector representing a ith voxel’s time series, i = 1,…,N. Since the mean of each voxel’s time series can be removed during regression, assume that each voxel has a mean of zero. Let the corresponding time series after global signal regression be represented by the column vector xi(t).
Regression of the global signal is accomplished by solving
where the global signal (a column vector) is given by
Estimation of the regression coefficient, βi, in Eq. (1) is given by
The average of the βi coefficient across all voxels is given by
Rearranging Eq. (4) provides the sum of the time series’ after global signal regression:
Inserting Eq. (3) gives
Let x1 be the seed voxel in a correlation analysis. Then the sum of dot products is given by
From Eq. (5) we can determine that
Then Eq. (6) becomes
This implies that the sum of dot products is less than or equal to zero:
The sum of correlation coefficients is given by
where μxi is the mean and σxi the standard deviation of the time series xi(t). Since we have assumed that μxi = 0 for i = 1,…,N, this equation reduces to
where M is the number of time points. If we assume that standard deviation values are unbiased across time series with negative and positive dot products, then the following equation holds:
Equation (7) demonstrates that if negative correlations don’t exist before global signal regression, they must be introduced by the technique. For all voxels that correlate positively with the seed, negatively correlated voxels must exist to balance the equation. This is compatible with the previous suggestion that global signal regression causes correlation strengths to be distributed around zero (Buckner et al., 2008; Vincent et al., 2006), although the equations do not provide any information about the resulting distribution. This theory section demonstrates that in fMRI data, after global signal regression, a seed voxel correlation analysis must find anti-correlations, regardless of the seed voxel chosen or the signal fluctuations it represents. This fact is a clear problem when trying to interpret the meaning of negative correlations after using global regression.
Simulations were performed to test the conclusions of the theory and to examine its practical implications. Two voxel and three voxel simulations were carried out using Matlab (MathWorks, Inc.). One thousand sets of time series were generated, each consisting of two and three time courses respectively. Resting state fluctuations were represented in each time series using a sine wave with a randomly chosen frequency between 0.1 Hz and 0.2 Hz and randomly chosen amplitude between 0 and 1. To this, Gaussian noise was added with zero mean and a standard deviation between 0 and 10. Each time series lasted 1000 time points with an effective TR of 1 second. For each set of time series, the global signal was regressed using the GLM method. In order to verify Eq. 3 and Eq. 5 empirically, the beta weights of the regression fits were averaged for each set and the sum of the time courses after global signal regression were calculated. The correlation coefficients between the first time course and the remaining one or two time courses in the set were calculated.
In another group of simulations, we investigated the effect of global signal regression on correlations between time series that exhibit similar resting state fluctuations that are slightly shifted in time. Null time series lasting 150 time points with a TR of 2 s were arranged on a 64×64 grid. A sine wave with a frequency of 0.1 Hz, amplitude of 1 and whose phase varied progressively along the x-axis from 0 to π/4 in equal increments was added to each time series to represent resting state fluctuations. Random Gaussian noise with a mean of 0 was added to each time series. The standard deviation of this noise was varied progressively along the y-axis from 0 to 10 in equal increments. The resulting dataset consisted of a 64×64 grid of time series with sine wave fluctuations of varying phase and increasingly large noise components. Global signal regression was performed by averaging over all 4096 time courses. Three correlation analyses were performed on the dataset, with and without global signal regression. The correlation of each time series with a sine wave of frequency 0.1 Hz and phase of zero was determined. Two 9×9 voxel regions-of-interest (ROIs) were drawn, one in a high and one in a low signal-to-noise ratio (SNR) region. The average time series in each ROI was determined and the correlation of the resulting time series with each voxel was calculated.
In the phase simulations above, every voxel contained a single frequency “resting state” fluctuation. A third “spatial extent” simulation investigated the influence of the fraction of voxels containing resting state fluctuations on correlation analyses when global signal regression is employed. One hundred artificial time series were generated lasting 300 s with a TR of 2 s. In each time series, physiological noise was simulated using a sine wave of 0.03 Hz (representing global signal changes, such as those induced by changes in depth of breathing) with random amplitude between 0 and 1 and random phase between 0 and π/4. Gaussian noise with a standard deviation of 1 was added. Resting state fluctuations were represented by a 0.1 Hz sine wave whose amplitude was 1 and 10 in low SNR and high SNR simulations respectively. This fluctuation was inserted in N of the 100 time courses, where N ranged from 1 to 99. For each iteration of the simulation, the global signal was calculated by averaging across the 100 time series and global signal regression was performed. A seed voxel time series was generated by regressing the same global signal from the 0.1 Hz sine wave resting state fluctuation. The average seed voxel correlation across the N voxels was determined along with the average across all remaining voxels. This whole process was repeated 100 times allowing study of the correlation values as a function of the percentage of voxels containing the fluctuations-of-interest.
This task combines both a breath hold block design and a visual block design task, varying the phase between the two in different conditions. The visual task produces a relatively small, focal activation in the visual cortex and was used as a model of resting state fluctuations. The global signal was modulated using a breath hold task that produces a large global BOLD signal increase (Kwong et al., 1992; Li et al., 1999). By combining the two, one can determine the effect of regression of a known global modulation on known focal “resting state” fluctuations (which are unknowns in actual resting state data). The extent to which the relative phase of the global and local signal changes influence local “resting state” connectivity measures was investigated.
Eight young healthy adults were scanned on a 3 T General Electric HDx MRI scanner equipped with a 16-element receive-only brain array coil under an IRB approved protocol after providing informed consent. Single shot, full k-space gradient recalled EPI was used for all functional scans. Whole brain coverage fMRI scans with 27 sagittal slices were collected with the following scan parameters: TR = 2 s, TE = 30 ms, matrix = 64×64, FOV/slice = 24 cm/4 mm, flip angle = 90°, reps = 150. Physiological data were recorded during each scan using a pulse-oximeter placed on the left index finger and a pneumatic belt positioned at the level of the diaphragm.
The subjects experienced five conditions: VisOnly, Synch, Synch + 10s, Asynch and RandVis (counterbalanced across subjects). In the VisOnly condition, an 8 Hz contrast reversing flashing checkerboard stimulus was presented using rear screen projection in a block design fashion 30 s OFF (grey fixation)/20 s ON. The subjects were asked to fixate on a cross in the centre of the screen. In the remaining conditions (Synch, Synch + 10s, Asynch and RandVis), the subjects were asked to perform a breath holding task. This consisted of a 30 s countdown presented centrally with paced breathing instructions, that is, “Breath In” and “Breath Out”, each lasting 2 s. When the countdown reached zero, the subjects are asked to hold their breath for 20 s with a similarly presented countdown. The task was timed to present “Breath Out” directly before the “Hold Breath” instruction so that breath holding occurred after expiration. In the Synch, Synch + 10s and Asynch conditions, a visual task identical to the VisOnly condition was presented together with the breath holding task with only the relative phase of the two tasks differing. In the Synch condition both the visual stimulus and breath hold onsets coincided temporally. In the Synch + 10s condition, the visual task was delayed by 10 s with respect to the breath holding task. Since responses to breath holds are typically delayed compared to BOLD activation responses, the aim of this condition was to cause the peak BOLD response for each task to overlap temporally. In the Asynch condition, the visual task was presented so that its ON period ended when the ON period for the breath holding task commenced. Finally, in the RandVis condition, an 8 Hz visual checkerboard, presented in an event-related design with a variable interstimulus interval, was displayed along with the breath hold task (30 s normal breathing, 20 s breath hold) such that for each second long interval there was a 50% probability of checkerboard presentation.
Data were pre-processed using AFNI (Cox, 1996). Reconstructed images were first corrected for motion using a rigid-body volume registration. Physiological noise correction was performed using the RETROICOR technique to remove unwanted signal oscillations at aliased cardiac and respiration frequencies (Glover et al., 2000). Data were time shifted to align separate slices to the same temporal origin. Voxel time series were bandpass filtered between 0.01 Hz and 0.1 Hz and converted into percent change by removing and then dividing by the mean. Global signal regression, using all voxels within the brain, was carried out using a GLM. Datasets with and without the global signal removed were carried forward into the next analysis step.
In each subject, a functionally defined visual cortex mask was generated using the VisOnly condition (without global signal removed) by correlating each voxel’s time series with the relevant task response and thresholding at P<1×10−6. This task-induced response was modelled by convolving the task timing with a Gamma-variate haemodynamic response function (irf(t) = k tr e−t/b where k is a scaling factor, r = 8.6 and b = 0.547; Cohen, 1997). Similar correlation analyses were performed on all conditions, with and without global signal removed, using the relevant visual task model. Correlation values were averaged across visual cortex and then across subjects. A seed region correlation analysis was performed by determining the 10 highest t-statistic voxels in visual cortex in the VisOnly condition (without global signal removed). The seed region time series was generated by averaging over these voxels in each condition. The correlation coefficient between this time series and the time course in each voxel of the brain was determined. For each subject, a histogram of these correlation values for each condition was calculated using voxels present in the visual cortex maps and also across all voxels in a whole brain mask. Group correlation maps were determined by converting each map to the standard stereotaxic coordinate system of Talairach and Tournoux (1988) using a 12-part linear transformation, averaging the arctanh of these maps and returning the tanh of the average (Fisher transformation).
Two resting state scans, each lasting 5 min, were collected on twelve young healthy adults with scanning parameters identical to the breath holding and visual data. Similar pre-processing steps were performed. After concatenation, the percentage change datasets, with and without global signal regression, were converted to the Talairach coordinate system as above. During this step, the data were resampled to 3 mm cubic voxels and spatially blurred using a Gaussian kernel with 6 mm full-width-at-half-maximum isotropic deviation.
Correlation maps were produced by extracting time courses from a seed region in the posterior cingulate/precuneus (PCC). This seed region was defined by drawing a 12 mm diameter sphere cantered around the previously published Talairach coordinate [−5, −49, 40] (Fox et al., 2005; Shulman et al., 1997). Time courses from voxels in this region were averaged and the correlation coefficient between this seed time course and each brain voxel was calculated. This process was carried out on the resting data with and without global signal removed. Group maps were obtained by averaging correlation values across subjects using a Fisher transform (for each voxel in Talairach space, the arctanh of the correlation value for each subject was averaged and the tanh of this average was returned). The individual subject correlation maps were thresholded at a correlation value of ±0.35 corresponding to a P-value of 9×10−10 (which, when corrected for multiple comparisons using a Bonferroni correction, equals P = 4.1×10−5). From the maps on which global signal regression was performed, the task-negative network (i.e. the network that deactivates during a task and is more active during rest) was defined on a subject-by-subject basis as those voxels that passed threshold and displayed a positive correlation with the PCC ROI. Since the PCC lies within the task-negative network, every region that positively correlates with it also lies within that network. On the other hand, the task-positive network (i.e. the network anti-correlated with the task-negative network) was defined as those voxels passing threshold displaying negative correlations with the PCC.
The correlation values in each of the task-positive and task-negative networks were compared in both the resting datasets with and without global signal removed. Since global signal regression is hypothesized to eliminate uninteresting low-frequency fluctuations from the data, the correlation values in the task-positive and task-negative networks were compared with a method that purports to remove low-frequency BOLD fluctuations due to changes in breathing depth. Respiration volume per time (RVT) correction was performed by determining the envelope of the respiration trace, subtracting the minimum from the maximum, dividing by the period of the respiration and then removing this time course (and time-shifted replicas) from each voxel using multiple linear regression (see Birn et al. for more details, Birn et al., 2006). Correlation values were averaged across the task-negative regions for each subject for the three conditions: without global signal or RVT correction, with RVT correction and with global signal regression. Similar averaging was performed in the task positive regions. Histograms of the correlation values across the entire brain, with and without global signal regression, were calculated.
Average time series in both the task-positive and task-negative networks were calculated on a subject-by-subject basis for datasets with and without regression of the global signal. Correlations between all combinations of these time series and the global signal were determined. The phases between these time series before global signal removal were calculated by determining the maximum correlation between both at multiple time lags. The percentages of voxels in the brain above threshold that positively and negatively correlated with the PCC region both before and after global signal removal were also calculated.
All 1000 sets of time courses in both the two voxel and the three voxel simulations demonstrate that Eq. (3) is true: the beta weights from the GLM average to exactly 1 (±~10−15) in all cases. Each time point summed across all time series in the set after regression of the global signal is exactly zero (±~10−13), demonstrating Eq. (5). The final correlation between the time courses after global signal removal in the two voxel case is exactly −1 (±~10−15) in all 1000 simulations. In the three voxel case, the sum of correlation values is always less than zero with a mean of −1.02 and standard deviation of −0.05 over the 1000 simulations. The maximum value of the sum was −0.83. This confirms Eq. (8) by demonstrating that after regression of the global signal, the sum of correlation values to a seed voxel is less than or equal to zero, regardless of the fluctuations it represents.
The phase simulation results are shown in Fig. 1. Correlations with the resting state sine wave demonstrate that the fluctuations are detectable in the majority of voxels even in low SNR time courses. Phase shifting the fluctuations by up to π/4 has little effect on detectability until global signal regression is performed, after which approximately half of the voxels display a negative correlation with the resting state sine wave. Correlations with averaged time series from ROIs show that the resting state fluctuations remain highly visible, whether the ROI is drawn from a high SNR (ROI1) or low SNR (ROI2) area, when global signal regression is not performed. After such a procedure, however, approximately half the voxels (52% and 48% for ROI1 and ROI2, respectively) display a negative correlation with the ROI. The phase simulations demonstrate that global signal regression can introduce negative correlations between time series with identical signal modulations that are slightly shifted in time.
The spatial extent simulations reveal that seed voxel correlation values are influenced by the size of the region containing the resting state fluctuations (see Fig. 2). Regressing the global signal from the data reduces correlation values as the spatial extent increases in both voxels containing fluctuations and voxels containing noise. Since noise voxels are not correlated with the seed voxel beforehand, global signal regression introduces negative correlations. In the low SNR case, this is a gradual decrease as spatial extent grows finally reaching a value of −0.3. In the high SNR case, this reduction is more rapid with the correlation values in noise voxels reaching this value when only 10% of voxels contain fluctuations. Both SNR cases are extremes and the situation in real data is most likely to be somewhere between the two. These simulations demonstrate that the reduction of correlation values due to global signal regression is dependent on the number of voxels containing fluctuations with greater spatial extent leading to a greater decrease. Not only does global signal regression reduce correlations in areas containing fluctuations, it introduces negative correlations in pure noise voxels.
The breath holding and visual data demonstrate that global signal regression can severely affect the outcome of resting state correlation analyses, particularly when the global signal (as represented by breath holding) and the “resting state” fluctuations (as represented by visual activation in this study) are highly similar. When BOLD signal increases due to visual activation and breath holding overlap temporally, as in the Synch + 10s condition (average correlation between visual task regressor and global signal was 0.61 ± 0.10), regression of the global signal severely reduces correlation values (see Fig. 3). In contrast, correlations are unaffected in the VisOnly condition since the global signal (which was not modulated my breath holding) is unrelated to the fluctuations-of-interest (correlation between task regressor and global signal: 0.08 ± 0.16). Similarly, when the global modulation (breath hold) is unrelated to the “resting state” fluctuations (visual activation), as in the RandVis condition, global signal regression only has a small impact on correlation values. The extent to which the global signal and the “resting state” fluctuations in the seed region are similar (i.e., the correlation value between the visual and global task) predicts almost exactly (r = 0.98) how much the global signal regression method will change the average correlation values in voxels containing the fluctuations.
Histograms of seed region correlation values across the visual cortex and across the whole brain are shown in Fig. 4. In all conditions, the distribution of correlation values in the visual cortex is skewed towards high values as one would expect. Global signal regression has little effect on the distribution of visual cortex correlation values in the VisOnly condition since the activation response and the global signal are unrelated. In the other conditions, global signal regression tends to widen the distribution towards lower values. Qualitatively, the greatest broadening of the distributions occurs in the Synch + 10s condition when the global signal is most correlated with the visual task regressor representing resting state fluctuations (0.61 ± 0.10), forcing correlation values to extend into the negative range. If global signal regression properly removes global manipulations, all distributions should resemble the VisOnly condition after performing the technique.
The whole brain distributions show that the breathing task causes a global modulation that skews the distributions when no global signal regression is performed. In all conditions, independent of the relationship between the global signal and the “resting state” fluctuations, global signal regression causes these skewed distributions to become more bell-shaped and centered on zero, thus introducing negative correlations. Assuming that global signal regression removes purely respiration related changes to reveal true neuronal fluctuations, the seed region connectivity maps after global correction should be identical, regardless of the breathing manipulation. Fig. 5 shows that large differences exist in the resulting correlation maps, none of which spatially resemble the VisOnly condition. Although voxels from the visual cortex are most positively correlated with the seed region in all conditions, areas that negatively correlate vary across breathing manipulations. The spatial organization of the negative voxels introduced by the correction is dependent upon the relative phase of the fluctuations-of-interest and the global signal. Coupled with the fact that highly correlated voxels remain highly correlated, this suggests that, after global signal regression, positive correlations may be more reliable than negative correlations. The extent to which these maps differ shows that global signal correction does not reliably remove the global artefact induced by the breath holding task.
Intrinsic correlations between a seed region in the PCC and all other voxels in the brain during a resting task are shown in Fig. 6. After global signal regression, approximately the same number of voxels are negatively correlated with the PCC seed region (task-positive network) as are positively correlated (task-negative network). When global signal regression is not performed, however, very few negatively correlated voxels pass threshold. The spatial distribution of the positive correlations is similar in both cases.
To probe the data more quantitatively, task-negative (correlated with the PCC) and task-positive (anti-correlated with the PCC) regions were defined for each subject. Average correlation values in these regions are shown for each subject in Fig. 7. Three conditions are depicted: the first draws values from data that have not been corrected for global fluctuations, the second from data that has undergone RVT correction and the third from data which has had the global signal regressed from it. RVT correction has been shown to remove some of the low-frequency fluctuations due to variations in depth of breath which causes spatially widespread confounds across grey matter and regions with large vessels including the PCC area (Birn et al., 2006). Unlike with global signal regression, correlations with the PCC did not greatly reduce when this correction was performed. Correlation values from the RVT corrected data are almost identical to those from the uncorrected data in both the task-negative and task-positive regions. On the other hand, global signal regression reduces correlation values in both areas. Reductions in the task-negative areas vary from 0.02 to 0.25 across all subjects. In task-positive areas, there is a greater variability across subjects in correlation values from the uncorrected data, ranging from −0.24 to 0.49 with a standard deviation of 0.23. The majority of subjects display strong correlations with the PCC in task-positive areas when global signal regression is not performed. After removal of the global signal, task-positive areas are strongly anti-correlated with the PCC and variability in the correlation values across subjects is greatly reduced (average correlation value is −0.42±0.07).
Histograms of the correlation values for each subject are shown in Fig. 8. The distribution is highly skewed towards positive values when global signal regression is not performed. Those subjects that show a negative correlation in task-positive areas in Fig. 7 correspond to the subjects with distribution peaks closest to zero. After global signal regression, the distributions revert to a more bell-shaped curve centred on zero. The tails of the distributions stretch further into the positive than into the negative correlation values. If correlation values in each voxel are compared on a voxel-by-voxel basis, before and after removal of the global signal, voxels on the high end of the distribution remain on the high end and voxels on the low end remain on the low end (data not shown). Qualitatively, comparing the distributions with those of Fig. 4, the resting state data appear to correspond most to the Synch + 10s condition. In this condition, the global signal was most highly correlated with the fluctuations-of-interest. If there are indeed two populations of voxels that are highly anti-correlated with each other, each including a large proportion of voxels in the brain, one would expect to see a bimodal-like distribution after global regression. A widened distribution in Subject 5 or the bump around −0.5 of Subject 9’s distribution may be evidence for this. However, the majority of subjects demonstrate no such behaviour suggesting the bell-shaped distributions of the correlation values are an artefact of the mathematics of the global signal regression procedure which forces negative correlations to appear.
Table 1 shows the sizes of regions that correlate both positively and negatively with the PCC before and after global signal regression. When the global signal is included, areas with a significantly positive correlation can constitute as little as 25% of voxels up to almost the entire brain (96%). Negatively correlated voxels are almost nonexistent with 0.59% being the greatest percentage in the twelve subjects. After global signal regression, positive correlations are reduced and negative correlations are increased (task-negative region = 9.65%±5.12%, task-positive region = 5.11%±4.65%). Phase differences between the task-negative and the task-positive time series before global correction range from 0.2 s up to 9 s. Comparisons with the global signal show that both the task-negative and task-positive time series correlate highly without removal of the global signal (0.75 ± 0.19 and 0.67 ± 0.18 respectively) and also display a large range of correlation values, from −0.51 to 0.62. After global signal regression, task-negative and task-positive time courses become highly negatively correlated with an average of −0.71 ± 0.14. The difference in correlation values before and after global signal regression spans a large range from 0.46 to 1.38.
The resting state data demonstrate that, when using the PCC as a seed region in a correlation analysis, anti-correlated networks are not significantly present until global signal regression is performed. As demonstrated in the breath holding and visual task, global signal regression renders the distribution of correlation values to be bell-shaped, centred on zero. After global signal regression, subject-to-subject variability in correlation values drawn from task-negative and task-positive areas is greatly reduced. Such regularized trends across subjects are a consequence of the global signal regression method.
The simple fact that any two time series become perfectly anti-correlated when their global signal is regressed from each should give cause for concern when interpreting anti-correlations. If the first time series contains a signal modulation of interest, the second will be negatively correlated with the first after global signal regression, regardless of the signal fluctuations it contains. Extending this concept to multiple time series, the theory demonstrates that global signal regression changes voxel time series such that approximately half become negatively correlated with a seed voxel. Due to the circularity of the technique (i.e., deriving a nuisance regressor from the data itself), all voxels could have similar interesting signal fluctuations but will display a range of seed voxel correlation values from highly positive to highly negative after its implementation. These findings alone should cause fMRI researchers to question the efficacy of this technique as a correction tool.
Resting state data are unique in that the researcher does not have a priori knowledge of task signal modulations. The simulations and breath holding data yield an insight into the effect global signal regression has on known fluctuations. For example, the spatial extent simulations show that time series in noise voxels outside an ROI containing resting state fluctuations become negatively correlated with those fluctuations when global signal regression is implemented. Not only does global signal regression reduce correlations in areas containing fluctuations, it introduces negative correlations in pure noise voxels. The degree of negative correlation depends on the extent of the ROI and the signal-to-noise ratio of the fluctuations. These simulations demonstrate that if there are large scale resting state networks in the brain, their influence on the global signal will force purely noise voxels to become negatively correlated with the resting network (e.g. the PCC region). When the SNR of the fluctuations is suitably high, only ~ 10% of the brains voxels need contain these signal modulations for the remaining voxels to become negatively correlated. The information in Table 1 suggests that in resting state data the task-negative network is large enough for this situation to occur.
The purpose of global signal regression is to remove any confounds that may mask neuronally-induced fluctuations. When confounds and resting state fluctuations are of similar frequency, as in the Breath Holding and Visual Data, their relative phases can severely affect the detectability of activation. In seed connectivity analyses where an explicit breathing confound exists, if global signal regression was completing the intended correction properly, one would expect the resulting connectivity measure distributions to resemble those of the VisOnly condition. However, this is not the case: global signal regression introduces negatively correlated voxels into the visual cortex. Even when no relationship between the global confounds and the fluctuations-of-interest exist (as in the RandVis condition), global signal regression does not sharpen the connectivity measure distributions as it should. These data demonstrate that global signal regression does not perform a correction that reveals the underlying neuronally-induced fluctuations.
In the Breath Holding and Visual Data, since the visual cortex represents such a small proportion of voxels it is expected that the seed voxels drawn from this area should bear little relation to most other voxels in the brain. When the breathing task introduces a global confound, whole brain connectivity measure distributions are skewed. Removing the global confound with regression results in distributions that are bell-shaped. This is unsurprising since the time series resulting from global signal regression are the residuals from a GLM analysis that assumes a Gaussian noise model. However, the spatial pattern of activations and deactivations varies depending on the breathing manipulation (see Fig. 5). This indicates that global signal regression does not remove the global breathing confound successfully. Although the visual cortex remains the most positively correlated area (as it should), negatively correlated areas can be situated in vastly different regions after global signal regression, suggesting that, in this case, the negative correlations are not neuronal in origin but an artefact of the global signal regression method.
The task-positive network is not evident in resting state data until the global signal is regressed from each time series. This is not a thresholding artefact; Fig. 7 and Fig. 8 demonstrate that lowering the threshold in the uncorrected map would not reveal negatively correlated voxels in the task-positive network. If it is possible to see positive correlations to the PCC before global signal regression, why are the negative correlations not visible? The standard answer is that the uninteresting, confounding fluctuations which are removed by the technique overshadow the underlying neuronally-induced fluctuations. However, removal of low-frequency confounds due to breathing depth changes using RVT correction barely alters correlation values. Even if RVT correction does not remove all global signals of interest, it is a global confound that is known to overlap well with the default mode network and correlates significantly with the global signal at CC = −0.5±0.13 (Birn et al., 2006). One would expect the RVT technique to somewhat reduce correlation values in the task-positive network.
There also exists a large amount of variability across subjects in the correlation values in the task-positive areas before any correction is performed. Comparing these values with the histograms in Fig. 8 shows that subjects with the lowest average correlation values in task-positive areas (see Subjects 4 and 11) have distributions that extend the furthest into negative values. Voxels that become the most strongly anti-correlated after global signal regression are the voxels that were least correlated prior to the regression. This suggests that when global signal regression distorts the distribution to become bell-shaped and centred around zero (as shown experimentally), voxels in the task-positive network always remain on the negative tail. Depending on the shape of this distribution, their average correlation value before correction can be both positive and negative. After global signal regression, the variability in these values is removed since all distributions are forced to be bell-shaped (thus satisfying the theoretical prediction that correlation values sum to less than zero). This implies that the drastic reductions of the correlation values in the task-positive network are due to the mathematics of converting a skewed distribution into a bell-shaped distribution centred on zero. Task-positive areas may be converted into negative correlations simply because they are the least correlated with the PCC to begin with.
If the global signal obscures fluctuations related to neuronal firing, why does the regression of this signal not uncover other positively correlated areas? Large areas of the brain are significantly positively correlated with the PCC before removal of the global signal. The task-negative network produced after global signal regression is always a subset of these positive areas. The areas that subsequently comprise the task-negative network can largely be determined by raising the threshold on the correlation maps without global signal regression (see Fig. 6). The global correction does not reveal any new areas that are positively correlated with the PCC. If the neurons firing in the PCC produced fluctuations that were not correlated to the global confound but were overshadowed by such a confound, we would not expect the positively correlated areas to be identical before and after global signal regression. This suggests that either the global signal is highly correlated to the PCC neuronal fluctuations or that fluctuations due to neuronal firing in the PCC are visible above the global signal confound. If the former were true then neuronal fluctuations in the PCC would be removed by global signal regression. If the latter were true, then anti-correlated networks should be visible above the global signal confound also. Since the task-negative network remains highly correlated after global signal regression and anti-correlated networks are not visible before, this suggests that a shifting of distributions by the technique is the most likely explanation for negatively correlated regions.
If negatively correlated voxels are an artefact of the correction technique, why are they in spatially contiguous areas and not randomly scattered throughout the brain? We have shown experimentally that global signal regression forces a skewed distribution to become bell-shaped, centred on zero and that voxels on the negative tail of this skewed distribution become the most negatively correlated. The reason that the task-positive network becomes anti-correlated after global signal regression is because it is the least correlated before. But why is it the least correlated? One possibility is that the correlation is reduced by a spatial heterogeneity in the “global” respiration related signal changes. Studies have shown that fMRI signal induced by variations in breathing depth (also reflected in variations in the end-tidal CO2) have spatial structure and are most prominent in areas associated with the default mode network (Birn et al., 2006; Wise et al., 2004). When seed voxels are chosen from the PCC, an area in the default mode network that is highly contaminated with breathing rate related confounds, areas outside the default mode network such as the task-positive regions will likely be the least correlated with the seed voxels. Also, phase differences can exist between time series drawn from the task-negative and task-positive areas and that these time series can have a low correlation value between them even though both are highly correlated with the global signal (see Table 1). This could be due to differences in vascular supply or haemodynamic properties that are unrelated to neuronal firing. The phase simulations show that this technique can introduce negative correlations between time series that display identical resting state fluctuations that are slightly shifted in time. The largest phase change in these simulations corresponds toa time delay of between 1.25 s and 12.5 s for frequencies between 0.01 Hz and 0.1 Hz. Although values on the upper end of this range are unlikely to be accounted for by haemodynamic delays alone, values on the lower end are within range of natural haemodynamic delay variation (Aguirre et al., 1998b; Handwerker et al., 2004; Saad et al., 2001) suggesting that areas whose neurons fluctuate simultaneously but whose haemodynamics vary temporally could display anti-correlations after global signal regression. Furthermore, underlying neuronal delays in low-frequency fluctuations on the order of seconds between networks would display high correlations with each other neuronally but would be deemed anti-correlated after global signal regression. (It should also be noted that pi/4 was an arbitrarily chosen phase delay and that if the maximum delay was smaller, similar correlation maps would be obtained). Spatial smoothing tends to increase the correlation of voxel time series with the global signal (Aguirre et al., 1997) and could also help consolidate sparsely positioned negatively correlated voxels into contiguous regions.
If global signal regression is a poor choice of correction method, which techniques should be used to remove noise confounds, be they physiologically related or not? The circularity of choosing a noise regressor from the data itself leads to the problems identified in this study. Confounding regressors derived from external measures, such as pulse oximeter, respiration belt and end-tidal CO2 readings, should not suffer from this problem. There exist other techniques that derive regressors from regions in the fMRI data, such as the ventricles, white matter, blood vessels and CSF, that may cause less of a problem. Studies have shown that anti-correlated networks do not exist after regression of movement, ventricle and white matter signals (Fox et al., 2008; Weissenbacher et al., 2008) and therefore these correction methods may not cause the same artefact as global signal regression. If a regressor is uncontaminated by the task in which we are interested (in this case the resting state fluctuations), then it may not introduce anti-correlations. However, the greater the correlation with the global signal, the larger a problem such a regression may cause by reducing true neuronally-related correlations and by increasing numbers of anti-correlated voxels. Using global signal regression to remove physiological noise clearly has confounding factors. While the optimal way to remove physiological noise isn’t clear, many of the best current methods require collection of cardiac and respiration time courses. Collecting such data should become a standard part of resting state connectivity studies.
This study has concentrated on correlation-based analyses, the most popular method employed to detect resting state networks. Other methods have been used to support the idea that the brain is organized into independent neuronal networks, for example, k-means clustering (Golland et al., 2008) or ReHo analysis (Long et al., 2008). However, removing the global signal is often included as a preprocessing step. An alternative technique that is now frequently used to study functional connectivity is independent component analysis or ICA (Birn et al., 2008; Calhoun et al., 2005; Kiviniemi et al., 2003; Kiviniemi et al., 2004; Long et al., 2008; McKeown et al., 2003; McKeown et al., 1998; McKeown and Sejnowski, 1998). With ICA, statistically independent maps and their associated time courses are extracted from the data. If the global signal is not removed before this technique is employed, it must be assigned to one or more components. Therefore, ICA may have an inherent global signal regression built into it. The time course for each component determines the behaviour of the network beyond all other components, including those that represent the global signal. Therefore, task-negative and task-positive networks determined by ICA may have associated time courses that artificially appear to be anti-correlated.
Global signal regression not only introduces anti-correlated regions in the brain, it also reduces positive correlation values. If such a simple mathematical process used as a pre-processing step affects correlation values so drastically, can we trust correlation techniques to give a true measure of connectivity in the brain? What do correlation values actually mean? In a clinical context, do small differences in connectivity measures between patient groups and controls really reflect deficits in neuronal connectivity? More importantly, can the global signal regression introduce false differences between these two groups? These questions must be taken into account before we can truly understand resting state neuronal networks in the brain, their relationships to each other and their importance in clinical conditions.
In this paper we have shown using theory, simulations, task and resting state data that anti-correlated networks in the brain may be a consequence of the global signal regression pre-processing step. Mathematically, this technique forces approximately half the voxels in the brain to become anti-correlated with a seed voxel. Global signal regression of data containing a known neuronal oscillation corrupted by a known respiration confound demonstrated that this technique is not successful in removing nuisance global effects and that the locations of anti-correlated areas are dependent on the relative phases of the global and seed voxel time series. In resting-state data, anti-correlated networks were not evident until global signal regression was performed. Other methods of low-frequency confound correction, such as RVT, have little effect on correlation values and, thus, do not reveal such networks. This study has demonstrated that employing global signal regression as a correction technique may cause spurious findings of negatively correlated regions in the brain. Great care should be exercised when interpreting such results.
This study was supported by the Intramural Research Program, National Institute of Mental Health, NIH. Thanks go to Phil Reiss and Clare Kelly for helpful discussions.