Topography and behavioral relevance of the global signal in the human brain

The global signal in resting-state functional MRI data is considered to be dominated by physiological noise and artifacts, yet a growing literature suggests that it also carries information about widespread neural activity. The biological relevance of the global signal remains poorly understood. Applying principal component analysis to a large neuroimaging dataset, we found that individual variation in global signal topography recapitulates well-established patterns of large-scale functional brain networks. Using canonical correlation analysis, we delineated relationships between individual differences in global signal topography and a battery of phenotypes. The first canonical variate of the global signal, resembling the frontoparietal control network, was significantly related to an axis of positive and negative life outcomes and psychological function. These results suggest that the global signal contains a rich source of information related to trait-level cognition and behavior. This work has significant implications for the contentious debate over artifact removal practices in neuroimaging.

In the imaging neurosciences, the global signal (GS) is defined as the timeseries of signal intensity averaged across all voxels in the brain, gray matter, or cortical gray matter. It is well known that non-neuronal sources including physiological noise caused by respiratory and cardiac events 1,2 and participant motion 3 contribute to the GS. As a consequence, GS regression became a pervasively adopted step in processing of resting-state fMRI data to attenuate these and other sources of noise 4,5 . However, in addition to containing artifactual information from various sources, the GS also contains information about ongoing neural activity 6 . Combined fMRI-electrophysiological studies in macaque monkeys permit analysis of spatiotemporal covariation between neural signal fluctuations measured with implanted electrodes and concurrent hemodynamic signals measured with fMRI. In one study, spontaneous fluctuations in local field potentials exhibited widespread positive correlations with fMRI blood oxygen level dependent (BOLD) changes over the entire macaque cortex 7 . More recently, neural origins of the global signal were indicated by inactivation of a neuromodulatory region of the basal forebrain, the nucleus basalis of Meynert. The nucleus basalis gives rise to the principal cholinergic as well as GABAergic projections to the cortex. Reversible pharmacological inactivation of the nucleus basalis in macaques resulted in regionally specific suppression of the global signal ipsilateral to the injection, further demonstrating a direct neuronal source of the global signal 8 .
Total estimates of baseline neuronal processing also come from magnetic resonance spectroscopy studies using 13 C radiotracers, which permit simultaneous measures of energy demand (CMR O2 ) in neurons and glia as well as neuronal activity as reflected by presynaptic release of the neurotransmitters glutamate and GABA 9 . Such work in rodents suggests that around 80% of neuronal energy in the cerebral cortex supports global neuronal activity at rest 10 . In light of these findings, Hyder and colleagues suggest that neither total baseline neuronal activity nor fluctuations in baseline neuronal activity can be neglected as merely representing non-neuronal factors 9 .
This contention is in line with positron emission tomography fluorodeoxyglucose studies in humans demonstrating that global signal amplitude is linked to changes in baseline glucose metabolism 11 . Additional evidence for neuronal contributions to the GS comes from EEG studies. The amplitude of the GS has been shown to relate to vigilance 12 and arousal 13,14 . Taken together, the electrophysiological, metabolic, and neuroimaging evidence clearly demonstrate that at least some neural information is carried in the GS.
A previous study demonstrated that GS correlations with each brain voxel follow a specific topography, with significant weights in the occipital lobe 3 . This spatial pattern is thought to reflect respiratory patterns, and has been exploited as a feature for artifact removal 15 . Other work has noted that GS correlations with each brain voxel are stronger in some functional brain networks than others, and they fluctuate in magnitude over time, as well as vary between individuals 16 . GS variability is higher in schizophrenia, and correlated with behavioral symptoms 17 . Differences in GS topographical representation have been documented between schizophrenia patients and controls 18 . Overall, while these previous findings hint at possible relationships between the GS and individual brain functional network architectures [16][17][18] , none of these studies have provided evidence of spatial structure in the inter-subject variability of global signal topography in healthy individuals.
A critical open question to address concerns how GS topography is related to human cognition and behavior. Here we systematically explored individual variation in GS topography and its relationship with individual differences on a range of cognitive and behavioral measures using canonical correlation analysis (CCA). We anticipated the existence of robust relationships between GS topography and phenotypic information, indicating functional relevance of this signal that is often discarded as noise.

Results
Descriptive Analyses: The GS manifests in a topographically specific manner. Capitalizing on the extensive neuroimaging and phenotypic de-identified data repository, the Human Connectome Project (HCP), we first assessed whether the GS has a distinctive spatial topography. Here, GS was operationalized by averaging the timeseries (t = 4800 time points from 4 scans for the majority of participants) across all surface-based cortical vertices for each participant. A GS beta map was then computed by regressing the GS at each vertex of each run of each subject (Fig. 1). The GS beta maps were then averaged across runs and across subjects yielding a mean GS Figure 1. Illustration of global signal beta map calculation. For each subject and each run, the global signal (X) was computed as the averaged timeseries across all cortical vertices. The global signal was then regressed from the timeseries of each vertex (Y), resulting in the GS beta map of a given run (β). Averaging across all runs within a given subject, we obtained the GS beta map of that subject. The GS beta maps for each subject were then used in subsequent PCA and CCA analyses (B).
testing for brain-behavior correspondence of the GS: canonical correlation analysis. To examine potential relationships between whole-brain GS patterns and subject-level behavioral and cognitive measures, we conducted a CCA between the 100 principal components derived from the GS topography data describing regional GS contributions for each subject, and 100 principal components derived from the HCP behavioral data. Non-parametric permutation analysis for null-hypothesis testing revealed a single statistically significant population correspondence between GS topography and behavioral profiles (r = 0.667, p < 0.001; Fig. 3). The pattern of canonical modes for the GS closely followed spatial patterns that recapitulate the frontoparietal control network (positive weights) and sensorimotor and visual networks (negative weights). In assessing the similarity between the first principal component (Fig. 2C) and the first canonical variate (Fig. 3B), a vertex-to-vertex correlation showed high convergence in topography (r = 0.81, p < 0.001). The pattern of canonical component weights for the cognitive and behavioral subject measures followed a general positive/negative gradient of life outcomes and psychological function, similar to that observed in a previous CCA of interregional resting-state functional connectivity and behavior in the HCP dataset 24 .
Robustness Analysis: Global signal topography and participant motion. Removal of the GS using linear regression has been proposed as the preferred method to attenuate participant motion related confounds in resting-state fMRI data 3,5 . Further, motion has been found to relate to individual differences in behavioral indices, suggesting that the propensity to move in the scanner (e.g. low task compliance), is a trait 25 . To confirm that the relationship between the GS and behavior we observed was not due to participant motion, we ran two more analyses: (1) computing the across-subject correlation between global signal beta estimates and mean head motion measures (DVARS, Framewise Displacement, FD) at each vertex, and (2) computing the across-subject correlation between canonical variate scores and mean head motion measures. As illustrated in Fig. 4, the brain topography associated with head motion is distinct from that observed to be associated with the positive/negative axis of behavior reported here (Fig. 3). In addition, no associations were observed between head motion (DVARS and FD) and the GS (r = −0.0002 and r = −0.0367, respectively) or behavior (r = −0.0001 and r = 0.0095, respectively). These supplementary analyses corroborate that meaningful associations between the GS and behavior are not likely to be mediated by head motion.

Aggressive Behavior
Antisocial Personality

Discussion
Resting-state fMRI has been widely embraced as a method to examine intrinsic functional brain networks 26 . Yet, researchers have acknowledged the challenges of separating neuronal and artifactual contributions in resting-state fMRI data 27 . The global signal (GS), or average signal intensity across the brain, is often removed from the timeseries via signal deconfounding to attenuate physiological and motion-related sources of noise 3,5 . However, the use of global normalization in fMRI data has long been debated 28 . Despite efforts towards consensus building with regards to the pre-processing step of removing the GS from fMRI timeseries via linear regression (GS regression) 29 , none has yet been attained 30,31 . Studies measuring electrophysiology and brain metabolism provide ample indication for some neuronal component in the GS [7][8][9] . In addition, fMRI work demonstrates that any regressors unrelated to true data noise, for example, also remove variance with network structure 32 . As such, implementing data processing procedures to remove the GS from a given fMRI timeseries may inadvertently discard relevant structured neural signal, particularly if the dimensionality of the data is low 33 . To date, there have been few systematic investigations into the relevance of the GS for studies attempting to link brain connectivity and individual differences in behavioral phenotypes 34 .
The present investigation examines individual differences in the spatial topography of the GS, followed by a CCA of these brain-wide GS patterns with the rich repertoire of behavioral data provided by the Human Connectome Project. First, we found that the subject-wise overall GS was related to cortical grey matter in a topographically-specific manner. The first principal component of the global signal, explaining ~6% of the variance, was reminiscent of the canonical frontoparietal control network 19,20 (positive weights) and sensorimotor and visual network (negative weights). The second component, explaining ~2% of the variance, represented the canonical anticorrelation between the default and dorsal attention networks (See Supplementary Fig. 1 for the first ten components of the GS).
Next, we conducted CCA using the principal components of the GS on the one hand, and phenotypic data on the other. This GS-behavior decomposition estimated using CCA aligns with a positive/negative axis of life experience and psychological indices. Greater frontoparietal control network weighting within the GS was significantly related to elevated scores on picture vocabulary, temporal discounting, life satisfaction and other measures on the positive axis, and lower aggressive antisocial behavior on the negative axis (Fig. 3). In fact, this positive/ negative axis of behavior is remarkably similar to that previously observed, with a critical difference that we used the subject-wise GS beta maps rather than inter-regional connectivity matrices of resting-state fMRI data in the CCA 24 . There appear to be many similarities between the current results and those of Smith and colleagues in the brain regions associated with this positive-negative axis of behavior. Key differences to note are brain functional connections robustly involving the cingulum and anterior temporal lobes 24 . In contrast, we found the strongest behavioral associations with lateral prefrontal cortical areas, anterior inferior parietal lobule, and other regions which comprise the frontoparietal control network 19 . The current results lead to the surprising observation that significant brain-behavior relationships can be derived from GS beta maps alone.
A possible explanation for our findings could be the underlying relationship between GS beta maps and functional connectivity. A vertex's value in the GS beta map represents how similar that vertex's timeseries is to the averaged timeseries from the whole brain. Mathematically, this GS beta value relates to the mean functional connectivity between the given vertex and all vertices, including itself. Given the results from Smith and colleagues, one might predict that the mean functional connectivity of each vertex is associated with the behavioral measures. However, it is still meaningful to show that individual differences in GS topography contain sufficient neural information to be associated with individual variation in behaviors. This suggests meaningful differences in functional MRI sources of variance (GS topography versus connectivity) that are associated with behavior.
The present study extends beyond the work of Smith and colleagues in multiple ways. First, we show similar brain-behavior associations with positive/negative axes as that reported by Smith and colleagues, although we have exclusively focused on the global signal aggregates that are commonly treated as a nuisance source. So, while that study examined what is commonly considered to be BOLD signal of interest, our analyses centered on what is often removed from the BOLD signal before performing any target analyses. We have determined how the global signal is differentially weighted across cortex (eg. global signal topography), with distinct magnitudes and variance. Most previous work that incorporates the global signal in the analysis uses one value per brain scan, which ignores the topographic specificity of this aggregate statistic. It is important to note that this regional weighting may impact the residual signal of interest when GS regression is applied.
The current findings contribute to two substantive research areas in contemporary network neuroscience, one theoretical and one methodological. First, our results speak to the theoretical systems neuroscience question of what aspects of neural activity are critically implicated in individual differences in cognition and behavior. Intrinsic functional connectivity MRI is a well validated approach to delineate and characterize functional-anatomic brain networks 35,36 . This method is also capable of revealing important dynamic aspects of neural processing that play an active role in cognition 37,38 . Intrinsic functional connectivity is thought to index stable individual features 39 , and has been associated with individual differences in a number of behavioral and cognitive domains as well as overall intelligence 40,41 . Such work has emphasized how interregional synchrony of brain activity involved in cognitive operations is associated with individual differences in performance of these cognitive operations 42 . The current work provides a novel complement to prior observations, suggesting that the GS topography alone carries structured information related to large scale brain networks and accounts for variability in behavior across individuals.
Methodologically, the current findings are relevant for the question of what effect specific pre-processing strategies have on the outcome of a given functional connectomics analysis. The practice of GS regression has come under scrutiny for a number of reasons. Even as early as PET and fMRI work from 1998, the validity of adjusting for effects of GS changes was called into question, as it may meaningfully alter results and thus interpretation of www.nature.com/scientificreports www.nature.com/scientificreports/ studies 28 . In resting-state fMRI analyses, the practice of GS regression can have additional unintended consequences. After GS regression, correlation values are mathematically centered on zero, which can produce spurious negative correlation values 43 . Researchers have noted that differences in caffeine intake can affect the GS such that caffeine leads to widespread decreases in connectivity and global signal amplitude, suggesting a neuronal source for GS that varies across individuals and time 44 . Another concern has been highlighted with simulated fMRI data, in which GS regression has been shown to artificially introduce correlations between brain regions and distort group differences in inter-regional correlations 33 . A study of autism spectrum disorder demonstrated that GS regression leads to a reversal in the direction of group correlation differences relative to other preprocessing approaches, with a higher incidence of both long-range and local connectivity differences that favor the ASD group 45 . Although GS regression has been shown to mitigate the effects of several sources of noise on estimates of functional connectivity 3,5,46 , the findings highlighted above suggest that functional connectivity values derived after GS regression has been performed may need to be interpreted with great caution.
One reason that deconfounding via GS regression persists despite these concerns is that the procedure minimizes the relationship between functional connectivity and motion 47 . Head motion is clearly related to some aspects of the GS, as we again demonstrate in the current work (Fig. 4). Data denoising and the process of dealing with motion artifacts is particularly challenging when working with large samples of "legacy" data of varying quality 35 . The practice of treating the GS as unwanted variation of no scientific value unfortunately remains in effect.
A recent protocol for mitigating head motion artifact in functional connectivity MRI suggests that GS regression is "singular in its ability to remove widespread artifact 5 ". Meanwhile, GS regression has been shown to strengthen the association of functional connectivity data with multiple behavioral phenotypes across cognition, personality and emotion 34 , possibly because the gain of having cleaner data outweighs the loss of neural information in the GS for some phenotypes. Therefore, the current findings demonstrating a relationship between GS topography and behavioral phenotypes are not contradictory to the strengthened associations reported by Li and colleagues. Depending on different research questions 29 , researchers should take great care in the choice of removing the GS during data preprocessing.
Some leaders in the area of fMRI data analysis have suggested that "…the field has reached a consensus that, as a pre-processing step, the global signal should not be removed 48 ". Many agree that it will be critical to continue to probe and better understand sources that contribute to the GS in future work 6,30 . There is growing evidence that prospective data acquisition procedures such as multi-echo fMRI can provide a means to effectively remove motion artifacts in fMRI data [49][50][51] . The current results do not speak to the question of whether or not GS regression should be utilized as a data processing strategy. Instead, we illustrate what tradeoffs can be anticipated when this step is included. Our results highlight the fact that a specific signal cannot always be unambiguously categorized as a confounder or not 52 . The nature of artifactual influences in the domain of imaging neuroscience may take a different form as a consequence of data richness, with hundreds of phenotypic variables often recorded for each individual. Our findings suggest that the neuroimaging community may have to carefully reconsider current deconfounding practices. Ultimately, the GS in human neuroimaging appears to be both signal and noise.

Methods
Dataset. Our population neuroscience study utilized the Human Connectome Project (HCP) S1200 release 53 .
Participants (N = 1094) were healthy young adults (ages 22-37) drawn from a population of twins and siblings. All imaging data were acquired on a customized Siemens 3 T Skyra at Washington University in St. Louis using a multi-band sequence. The structural images were 0.7 mm isotropic. The resting-state fMRI data were 2 mm isotropic with TR = 0.72 s. Two sessions of rs-fMRI data were collected on consecutive days for each subject, and each session consisted of one or two runs. The length of each rs-fMRI scan was 14.4 min (1200 frames). Details of the data collection can be found elsewhere 53,54 . Informed consent was obtained from all subjects. Details about behavioral measures can be found in HCP S1200 Data Dictionary and 55 . All methods were carried out in accordance with relevant guidelines and the University of Miami Institutional Review Board approved the study.
Motion censoring was performed to remove subjects and runs with high motion. The to-be-censored frames (outliers) were identified by three stages. First, volumes with FD > 0.2 mm or DVARS > 75 were marked as outliers. Second, one frame before and two frames after these volumes were also flagged as outliers. Finally, remaining segments of data that lasted fewer than five contiguous volumes were censored. All censored frames marked in these three steps were discarded when computing global signal beta maps. BOLD runs with more than half of the frames flagged as censored were removed, resulting in 1030 subjects after motion censoring. Note that the FD threshold we picked is more conservative than that used in previous literature 11,46,60,61 . The DVARS was selected to achieve a similar number of censored frames as flagged by the FD threshold. In addition, two subjects did not have adequate family structure information (as used in the permutation testing procedure described below), leaving 1028 subjects for the final CCA.
Delineation of global signal beta maps. The procedure to generate global signal beta maps is shown in Fig. 1. For each subject and each run, the global signal was calculated by averaging the timeseries across all cortical vertices in grey matter. Then for each run, the global signal and a vector of ones were regressed together (i.e. with bias term/intercept) from the fMRI signal for each vertex in the fs_LR space using ordinary least squares. To exclude the effects of high motion, censored frames were ignored when univariate linear regression was performed to regress the overall GS against each voxel's BOLD activity fluctuations 34,62 . The ensuing regression (beta) coefficient corresponds to the differences in the global signal for each vertex. This constituted the global signal beta map for each scan. Since global signal beta maps were computed from multiple scans for each subject, we averaged the maps from multiple scans for each subject. The averaged beta map captures the overall relationship between the global signal and the BOLD signal at each vertex of the examined subject. The average global signal beta map for each subject was used for further analyses.
It is worth noting that there are multiple ways of computing the GS in the literature. While we computed the GS using just the cortex 34,63 , one can also compute the GS using the whole brain, or all gray matter voxels 51 . However, all definitions of GS lead to very similar results. For example, in the current data, the average correlation between GS computed from all gray matter locations (i.e., grayordinates) and the GS computed from just the cortical surface was 0.95. Thus, it is unlikely that alternate definitions of the GS would significantly alter our results. canonical correlation analysis. To interrogate the relationship between the GS and behavior, we conducted a CCA. CCA is a natural choice of method because this machine-learning algorithm computes linear combinations of the original variables for each of two multivariate datasets that, together, maximize the linear correspondence between both variable sets. Each ensuing canonical mode, reflecting a pair of canonical variates, is indicated by a linear weighting of behavioral measures and a linear weighting of cortical vertices that are maximally correlated with each other. The strength of association (r) between the two variable sets was tested for robustness and statistical significance using a statistical null-hypothesis testing permutation framework 24 . Behavioral data exclusion and pre-processing. Several behavioral measures contain large amounts of missing observations, have heavily skewed distributions, and/or may represent potential confounds. To ensure these variables did not adversely affect our analysis, we followed a similar data-exclusion and pre-processing procedure as that described in 24  In total, 143 behavioral measures were excluded (including confounds). Following removal of these variables, all categorical variables were dummy-coded, and all variables were then z-score normalized. This resulted in a battery of 177 standardized behavioral measures for the CCA analysis. See Supplementary Fig. 2 for analyses without behavioral data exclusion.

Data reduction.
To avoid overfitting in the CCA analysis, and to contend with the high dimensionality of the datasets, we conducted a principal component analysis (PCA) individually on the pre-processed behavioral measures and on the global signal beta maps (z-score normalized), analogous to previous work 24 . Since we obtained a global signal beta map for each subject, we treated every subject as an observation and each vertex as a feature dimension in PCA. For example, the first principal component of global signal beta maps captured the most inter-subject variance in the global signal spatial topography. 100 principal components were derived from each of the behavioral and global signal beta maps. As in 24 , the covariance matrix of the behavioral measures was projected onto the nearest positive-definite covariance matrix, which avoided any imputation of missing observations.
A CCA was then conducted on the 100 behavioral and 100 global signal principal component scores. CCA estimates canonical variate pairs, or linear combinations of the behavioral with its links to a linear combination of global signal principal components, that are maximally correlated. Canonical variate pairs are organized in terms of decreasing magnitude, meaning the first canonical variate pair explains the largest amount of variation in the data, followed by the mode with second-largest explained variance, and so forth. The maximal number of canonical components (N = 100) was estimated. The CCA output of interest was the overall strength of the correlation between each canonical variate pair; the pattern of canonical weights for the original behavioral measures, on the one hand, and cortical vertices on their canonical variate, on the other hand. The canonical weights for each measure and cortical vertex were computed as the correlation between the canonical variate scores and the original behavioral measures and GS beta values across subjects.

Statistical significance testing of canonical variate correlations (r).
To test for the statistical significance of the correlation strength (r) for all discovered brain-behavior associations (e.g. canonical variate pairs), we carried out the null permutation testing framework used in 24 . This flexible testing scheme aimed at rejecting the null hypothesis that the obtained pairs of canonical variates are likely due to random noise linking our brain and behavior data. The permutation test proceeded as follows: (1) subjects (i.e. rows) behavioral principal component scores were randomly re-ordered, respecting family structure 64 , to intentionally break the dependence structure between the two variable sets, (2) the CCA analysis was re-run between the original GS principal component scores and the re-ordered behavioral principal component scores, and (3) the first canonical variate pair's correlation coefficient (the maximum possible correlation coefficient) was placed into a null distribution of correlation coefficient values supposing a lack of brain-behavior association, (4) this process was repeated 10,000 times, and (5) the original canonical variate pairs were declared statistically significant if their associated correlation coefficients exceeded the 99.9% percentile (i.e., p < 0.001 of the null distribution defined above). As noted in Results, only the first canonical variate pair was statistically significant (p < 0.001).

Data Availability
The MRI and behavioral datasets used in this study are available to the public from the Human Connectome Project (S1200 release; https://www.humanconnectome.org/study/hcp-young-adult/document/1200-subjects-data-release).

code Availability
All code used to analyze data and create figures for the current study are available with no restrictions. Customized Matlab code can be accessed here: https://github.com/tsb46/Global-Signal-CCA-Analysis.