Abstract
A major focus of brain research recently has been to map the resting-state functional connectivity (rsFC) network architecture of the normal brain and pathology through functional magnetic resonance imaging. However, the phenomenon of anticorrelations in resting-state signals between different brain regions has not been adequately examined. The preponderance of studies on resting-state fMRI (rsFMRI) have either ignored anticorrelations in rsFC networks or adopted methods in data analysis, which have rendered anticorrelations in rsFC networks uninterpretable. The few studies that have examined anticorrelations in rsFC networks using conventional methods have found anticorrelations to be weak in strength and not very reproducible across subjects. Anticorrelations in rsFC network architecture could reflect mechanisms that subserve a number of important brain processes. In this preliminary study, we examined the properties of anticorrelated rsFC networks by systematically focusing on negative cross-correlation coefficients (CCs) among rsFMRI voxel time series across the brain with graph theory-based network analysis. A number of methods were implemented to enhance the neuronal specificity of resting-state functional connections that yield negative CCs, although at the cost of decreased sensitivity. Hubs of anticorrelation were seen in a number of cortical and subcortical brain regions. Examination of the anticorrelation maps of these hubs indicated that negative CCs in rsFC network architecture highlight a number of regulatory interactions between brain networks and regions, including reciprocal modulations, suppression, inhibition, and neurofeedback.
Introduction
Over the last two decades, resting-state functional connectivity (rsFC) assessed with functional magnetic resonance imaging (rsFMRI) has become an important tool in the study of brain function and pathology. Interest in rsFC networks was catalyzed by initial findings of anticorrelations between the default mode network (DMN) and the task-positive network (TPN) (Fox et al., 2005), which mimicked the reciprocal modulation of brain activity between these two networks during cognitive and attention-demanding tasks (Gusnard and Raichle, 2001). However, anticorrelations between DMN and TPN, found by Fox et al. (2005), were later established to be artifacts arising from the employment of global signal regression (GSR) in the data analysis pipeline (Murphy et al., 2009). A few studies (Chang and Glover, 2009; Fransson, 2005) have reported the persistence of anticorrelations between DMN and TPN (Chang and Glover, 2009; Fransson, 2005) and between DMN and the salience network (Uddin et al., 2009) even when GSR was not employed during data analysis. However, in general, negative cross-correlation coefficients (CCs) between different brain regions have been found to be weaker and less reproducible compared with positive CCs (Chang and Glover, 2010; Shehzad et al., 2009) when GSR is not performed.
Some studies have reported that negative CCs occur between regions with low direct anatomical connectivity (Honey et al., 2009). However, recent reviews (Buckner et al., 2013) have shown that rsFC can occur between regions with no direct anatomical connections. Furthermore, rsFC is considered to represent both anatomically constrained (correlations between regions, which have direct anatomical connections) and state-dependent signal couplings (Buckner et al., 2013).
Anticorrelation between brain regions can reflect a number of state-dependent regulatory interactions. One example of these is reciprocal modulations in the activity between two separate brain functional networks (e.g., DMN and TPN or cognition and emotion). These modulations are thought to be mediated through indirect anatomical connections between the two networks (Buckner et al., 2013; Drevets and Raichle, 1998; Gusnard and Raichle, 2001). Another example is active suppression of a brain network in the resting state through anticorrelations between specific components of the network mediated by direct anatomical connections. The thalamus and basal ganglia, which have a regulatory role in all cortical brain networks (Crosson, 2013; Postuma and Dagher, 2006), are examples of structures that could be hubs of anticorrelation (HACs) resulting from regulatory mechanisms.
One reason for the observed weaker negative CCs compared with positive CCs between brain regions could be that the latter are artificially inflated due to common noise processes across the brain, for example, residual or unmodeled physiological noise (Bianciardi et al., 2009) and global signal (Scholvinck et al., 2010). Anticorrelations between two networks (e.g., DMN and TPN) can also be diminished by dynamic changes in the brain states during a typical rsFMRI scan (Chang and Glover, 2010). Low reproducibility of anticorrelated networks across a group could also reflect variability in the location of drivers of anticorrelations in rsFC. In general, one can surmise that anticorrelations are more likely to reflect state-dependent signal couplings, whereas positive CCs are more likely to highlight anatomically constrained spontaneous fluctuations. Consequently, it is advisable to examine rsFC networks derived from negative and positive CCs separately.
In this preliminary study, we examined the hubs of rsFC characterized by high degree centrality (DC) in two separate whole-brain high-resolution rsFC graphs formed with rsFC connections determined by positive and negative CCs for each subject. The data analysis pipeline included a number of steps to enhance the neuronal specificity of the rsFMRI signal as well as to increase reproducibility of anticorrelation maps across subjects. We did not employ any intensity-based nuisance variable regression (e.g., global signal), which could lead to artifactual anticorrelations (Murphy et al., 2009). We also restricted the signal to gray matter and masked out voxels near high susceptibility contrast boundaries. Furthermore, to overcome reductions in reproducibility of negative CC maps across subjects (Shehzad et al., 2009) due to anatomical variations, the location of drivers of anticorrelations in a given brain area (e.g., supplementary motor area [SMA]) were allowed to vary by 5 mm across subjects. Our hypothesis was that the HACs based on negative CCs and their associated rsFC anticorrelation maps will shed light on different regulatory interactions between networks and areas that subserve brain processing during the resting state.
Materials and Methods
Participants
Twelve healthy right-handed participants were recruited from the Atlanta metropolitan area (7 male; median age=27 years; native language=English). Exclusion criteria for the study included pregnancy, metallic implants, or other contraindications to MRI; history of psychiatric or neurologic illness or substance abuse. Written informed consent was obtained from participants under an IRB-approved protocol.
Image acquisition
MRI images were acquired on a 3T Siemens TIM Trio scanner with a 32-channel receiver array head coil. Participants were instructed to keep their eyes open and blink at a normal rate during the resting-state scans (Gopinath et al., 2011). rsFMRI scans were acquired with a multiband slice-accelerated gradient-echo EPI sequence (Moeller et al., 2010) with imaging parameters as follows: field of view (FOV)=220 mm; repetition time (TR)/echo time (TE)=2000 / 24 msec; flip angle (FA)=45°; phase partial Fourier factor=5/8; 48 interleaved 3×3 mm resolution axial slices of 3 mm width covering the whole brain; and 270 scan volumes. The relatively short TE, partial Fourier fraction, and thin slices were prescribed to reduce magnetic susceptibility artifacts. In addition, a whole-brain three-dimensional (3D) T1-weighted MPRAGE sequence (FOV=240 mm; TR/TI/TE/FA=2250/900/3 msec/9°; 0.9×0.9×1 mm resolution) was acquired to provide anatomical reference. Foam padding was provided to minimize head motion. A respiration belt and finger pulse plethysmograph (Biopac, Goleta, CA) were used to acquire physiological waveforms time-locked to fMRI data acquisition.
Data analysis—preprocessing and physiological noise reduction
Data analysis was performed using the AFNI and FSL software packages, in-house programs written in Matlab (Natick, MA), and the Brain Connectivity Matlab toolbox (Rubinov and Sporns, 2010). Voxel time series in the rsFMRI dataset temporally shifted to account for differences in slice acquisition times. The 3D scan volumes were then registered to a base volume to account for motion. The resultant rsfMRI time series were detrended of the signal related to cardiac and respiratory oscillations at their respective main frequencies, first harmonics and low-frequency fluctuations in BOLD signal arising from changes in respiration volume per time and in the heart rate, using the most conservative physiological noise removal technique (Bianciardi et al., 2009). The resultant time series were coregistered to the T1-weighted high-resolution anatomical scan and spatially normalized to the MNI152 template. The spatially normalized time series data were low-pass filtered (cutoff frequency=0.1 Hz) to better isolate low-frequency BOLD rsFC fluctuations. To enhance neuronal specificity of the rsFMRI data, a subject-specific brain mask was created. The brain mask was restricted to voxels, which contained at least 80% gray matter (segmented using the high-resolution anatomic). Furthermore, voxels in this mask at brain edges and high susceptibility contrast boundaries were removed from analysis. The rsFMRI data for each subject were spatially smoothed with an isotropic Gaussian filter (FWHM=5 mm) restricted to its subject-specific brain mask.
Graph-based hubness analyses
The preprocessed rsFMRI voxel time series data were parcellated into nodes to form rsFC network graphs. We adopted a high-resolution node parcellation approach. Each voxel in the subject-specific brain mask acted as a node. For each subject, separate rsFC network graphs (with identical nodes) for positive and negative CCs were constructed to overcome the inherent bias in positive CC distributions. The negative CC graph employed a binary connectivity distance matrix thresholded at CC<−0.12, which corresponds to an uncorrected connection level, p<0.05, to ensure rejection of spurious connections. The positive CC graph employed a binary connectivity distance matrix whose threshold CC was set for each subject individually by ensuring sparsity=0.3 for each subject in the group. This threshold is roughly the minimum sparsity at which small-world properties are maintained in the graph (van den Heuvel et al., 2008). For each node, the DC was calculated and hubness was determined by standardizing the DC estimate across each subject's dataset: Hi=(DCi – mean(DC))/std(DC), where Hi is the standardized hubness of node i (Buckner et al., 2009). To account for the hypothesized variability in the location of HACs across subjects, group-level hubness maps were constructed by first replacing the hubness score of a given voxel in each individual subject hubness map by the maximum of hubness scores of voxels within a 5 mm distance to it in the MNI space. Hubness statistical parametric maps (for negative and positive CC graphs) were constructed through one-sample t-tests on the respective individual subject hubness maps.
Finally, to further probe the functional connections characterized by anticorrelations, rsFC maps were constructed for carefully selected HACs (see the Results section for description), through one-sample t-tests on z-transformed CC (with CC<0) maps of the corresponding nodes. To account for variability in the location of the drivers of anticorrelations as well as variability in the location of gray matter gyri across subjects, the z-transformed CC of each voxel was set as the maximal negative z-statistic within a 3 mm distance of it in the MNI space. The group t-test maps generated from each of the above analyses were clustered and the significance of the results accounting for multiple comparisons were derived by means of Monte Carlo simulation of the process of image generation, spatial correlation, intensity thresholding, masking, and cluster identification (Forman et al., 1995). The significance of hubness and rsFC maps reported in the Results section are corrected for multiple comparisons by this procedure.
Results
Distribution of CCs in subject-level rsFC network architecture
None of the subjects examined displayed more than 2 mm motion during their rsFMRI scan, barring a few isolated (less than 1%) volumes, which were censored during analysis.
Figure 1 shows distribution of CCs for two representative subjects. An average of 20% (standard deviation=10%) of rsFMRI connections in the brain connections were characterized by negative CCs. The negative side of the CC distributions was not altered significantly (at p=0.05; considering the mean and median of the distribution) after removal of physiological noise from the voxel time series. On the other hand, the positive CC distribution was significantly (p<0.001) shifted to the left (i.e., CCs were reduced) after physiological noise removal. The average sparsity in negative CC graphs across subjects for a CC threshold = −0.12 (uncorrected p<0.05) was 0.09. For the positive CC graphs, the average CC corresponding to selected sparsity threshold (0.3) was CC=0.30. Sparsity of 0.1 in positive CC graphs corresponded to average CC threshold=0.48.

Distribution of cross-correlation coefficients (CCs) of the correlation matrix for two representative subjects.
Location of HACs in subject-level rsFC network architecture
Table 1 lists the hubness z-scores of HACs in three representative regions: SMA, thalamus, and insula. The location of HACs in these and other observed regions, for example, putamen, paracentral lobule (PCL), and subgenual cingulate (SgC), varied by around 5 mm across the group. S10 did not exhibit HACs in any region.
Location of Hubs of Anticorrelation (Hubness Maxima) in Selected Regions in Individual Subjects
SMA, supplementary motor area; VL thalamus, ventrolateral thalamus.
Group-level HACs
Significant (p<0.05) HACs (see Fig. 2 and Table 2) were found bilaterally in a number of regions, including the insula and superior temporal cortex (STC), SMA, PCL, SgC, and ventrolateral (VL) prefrontal cortex. Furthermore, the thalamus, basal ganglia, hippocampus, and amygdala were densely populated with HACs. Large parts of the thalamus (bilateral: VL, ventral anterior [VA], ventroposterio medial, medial dorsal [MD], and pulvinar), putamen (bilateral caudal and left rostral) as well as hippocampus acted as HACs.

Group average hubness (H) z-score map for significant hubs of anticorrelation (HACs) in brain resting-state functional connectivity (rsFC) network architecture. These hubs were found to be significant (cluster level p<0.05) across the group in a one-sample t-test on H-maps. Slice locations in MNI coordinates. The seed locations for the rsFC maps shown in Figures 3 –6 are indicated with arrows: supplementary motor area (SMA); paracentral lobule (PCL); subgenual cingulate (SgC); ventrolateral thalamus and caudal putamen.
Location of Hubs of Anticorrelation in Significant Clusters in the Group-Level Hubness t-Statistic Map (Multiple Comparisons Corrected Cluster Level p<0.05)
FEF, frontal eye fields; VPM/VPL thalamus, ventral posterior medial, or ventral posterior lateral thalamus.
Group rsFC maps of selected HACs
Based on the location of HACs described above, we proceeded to examine the rsFC maps based on negative CC of selected HACs in the cortex, basal ganglia, and thalamus to further probe the mechanisms that are highlighted by the anticorrelations in rsFC network architecture. HACs in SMA and PCL were chosen since they are part of well-known motor programming and somatosensory perception and pain (DeLong and Wichmann, 2009; Favilla et al., 2014). SgC was chosen since it was part of the emotion and limbic networks implicated in major depressive disorders (Johansen-Berg et al., 2008). HACs in basal ganglia and thalamus were examined since these regions are involved in all brain networks (Postuma and Dagher, 2006). Specific locations in the caudal putamen and VL thalamus were chosen based on the strength of their hubness score (maxima) in basal ganglia and the thalamus, respectively.
The HAC in left SMA (Fig. 3) showed significant (p<0.05) anticorrelations bilaterally with the caudal and rostral putamen, ventral striatum, thalamus, and cerebellum and in the left hemisphere with the insula, amygdala, and hippocampus head. Right PCL (Fig. 4) exhibited anticorrelations bilaterally with the dorsal anterior cingulate, dorsal striatum, caudal putamen, thalamus, and cerebellum and with the left hemisphere insula and STC. Left SgC (Fig. 5) exhibited anticorrelations bilaterally with the MD thalamus, pulvinar, caudal putamen, SMA, PCL, insula, and STC and the left hemisphere amygdala, hippocampus head, and ventral striatum. Subcortical hubs in the right caudal putamen and right VL thalamus exhibited (Fig. 6) anticorrelations with all parts of the cortex. The right caudal putamen also exhibited anticorrelations bilaterally with the thalamus, hippocampus, and caudate, and left caudal putamen and amygdala, and right rostral putamen. The right VL thalamus exhibited anticorrelations bilaterally with all areas in basal ganglia, hippocampus, amygdala, and thalamus, except the VL thalamus nuclei.

Group t-test (rsFC) map of brain regions exhibiting significant (cluster level p<0.05) anticorrelated functional connections to the strongest HAC in supplementary motor area. Slice locations in MNI co-ordinates.

Group t-test (rsFC) map of brain regions exhibiting significant (cluster level p<0.05) anticorrelated functional connections to the strongest HAC in the paracentral lobule. Slice locations in MNI coordinates.

Group t-test (rsFC) map of brain regions exhibiting significant (cluster level p<0.05) anticorrelated functional connections to the strongest HAC in subgenual cingulate cortex. Slice locations in MNI coordinates.

Group t-test (rsFC) map of brain regions exhibiting significant (cluster level p<0.05) anticorrelated functional connections to the strongest HAC in
Group-level hubs of positive correlation
Significant (p<0.05) hubs based on positive CCs were found (see Table 3) in a number of bilateral areas, including parts of the DMN (precuneus/posterior cingulate, retrosplenial cortex, lateral parietal cortex, and inferior temporal region), cerebellum, occipital, and somatosensory and motor cortices, consistent with previous studies (Buckner et al., 2009; Cole et al., 2010). The hubness z-scores were comparable with the values reported in (Buckner et al., 2009).
Location of Hubs Based on Positive Correlation in Significant Clusters in the Group-Level Hubness t-Statistic Map (Multiple Comparisons Corrected Cluster Level p<0.05)
Discussion
HACs and mechanisms subserving associated anticorrelation maps
Examination of anticorrelation maps of individual cortical HACs (SMA, PCL, SgC) in task-positive and emotion networks revealed (Figs. 3 –5) that resting-state fluctuations between different constituent areas of each of the sensorimotor, pain, and emotion networks (Crosson, 2013; DeLong and Wichmann, 2009; Drevets and Raichle, 1998; Favilla et al., 2014; Postuma and Dagher, 2006) are anticorrelated with respect to each other at rest. The functional connections highlighted by these anticorrelations are subserved by well-known direct anatomical pathways connecting these networks (Johansen-Berg et al., 2008; Vergani et al., 2014). The motor initiation network consists of the primary motor cortex (M1), SMA, insula, cerebellum, caudal putamen, and VL thalamus (DeLong and Wichmann, 2009). The RsFMRI signal from parts of SMA shows anticorrelations between parts of all these areas (Fig. 3). The phenomenon of anticorrelated rsFMRI fluctuations between components of the motor initiation network may support active suppression of this network at rest. Specifically, the anticorrelations may reflect the case where, when the parts of SMA around the HAC become active, other parts of the motor initiation network are suppressed to prevent motion during the resting-state scan. Similarly, the HAC in PCL exhibited (Fig. 4) significant anticorrelation with parts of the anterior cingulate, anterior insula, putamen, thalamus, amygdala and medial temporal lobe, and cerebellum at rest. These areas are part of somatosensory and pain networks (Favilla et al., 2014; Postuma and Dagher, 2006) and the observed anticorrelations may serve to suppress excitability of these networks during the resting-state scan. Similarly, the HAC in SgC (Fig. 5) exhibited anticorrelations with parts of the emotion network (Drevets and Raichle, 1998; Johansen-Berg et al., 2008)—amygdala, ventral striatum, MD thalamus, hippocampus, and anterior insula—perhaps reflecting processes subserving the suppression of excitability of this network during the neutral affective condition engendered by the rsFMRI scan.
Interestingly, the rsFC (Fig. 5) map of the HAC in SgC also showed significant anticorrelations with SMA, PCL, posterior insula, and caudal putamen—all components of motor and somatosensory networks. Analogously, the HAC in SMA (Fig. 3) exhibited anticorrelations with the ventral striatum, amygdala, and hippocampus head—all components of the emotion and limbic networks. Thus, some of the anticorrelations in the rsFC network architecture may also reflect reciprocal modulations between sensorimotor and emotion networks (Drevets and Raichle, 1998; Gopinath et al., 2011).
Dense arrangements of strong HACs were found in the thalamus and basal ganglia (Fig. 2). Anatomical projections exist from basal ganglia and thalamus to all parts of the cortex (Johansen-Berg et al., 2005; Tziortzi et al., 2013), and basal ganglia–thalamocortical circuits are involved in all brain functions (Crosson, 2013; Postuma and Dagher, 2006). Thus, it is not surprising that HACs in the thalamus and putamen exhibit anticorrelation (Fig. 6) with all areas of the cortex. This is also consistent with the presence of the thalamus and basal ganglia in anticorrelation maps of cortical HACs (Figs. 3 –5). Interestingly, the HAC located around the right VL thalamus nucleus exhibited anticorrelations with parts of other thalamic nuclei (e.g., bilateral VA, MD, and pulvinar nuclei). Since thalamic nuclei do not possess monosynaptic connections to each other (Crosson, 2013), anticorrelations between thalamic nuclei must be caused by polysynaptic connections through the cortex, that is, through interactions between different functional networks. These interactions could be a reflection of regulatory mechanisms between different functional networks (e.g., suppression/inhibition, reciprocal modulations, or neurofeedback). The VL thalamus is a motor nucleus, whereas VA-, MD-, and pulvinar thalamus are primarily involved with cognitive processing (Crosson, 2013; DeLong and Wichmann, 2009), with MD thalamus also connected to the limbic system (Postuma and Dagher, 2006). Thus, anticorrelations between thalamic nuclei may indicate that if one thalamic nucleus is activated, the functions associated with other thalamic nuclei will become less active. All the above hypothesized mechanisms for anticorrelations in rsFC network architecture deserve further attention as they may indicate how different network components will interact under different situations.
Distribution of CCs and effects of vasomotion
At an average, about 20% of rsFC connections in the brain were characterized by negative CCs. The conservative GM masking in conjunction with detrending of signal related to respiration and cardiac waveforms minimized the CSF, vasomotion, and physiological contributions to the signal. The careful removal of voxels at the brain edges and high susceptibility contrast boundaries also enhanced the neuronal specificity of the rsFC networks. Negative CC graphs were very sparse compared with positive CC graphs. At the negative CC threshold chosen, of CC=−0.12 (uncorrected p<0.05), the average sparsity across all datasets was 0.08. On the other hand, sparsity of 0.1 corresponded to a relatively high average CC threshold (CC≥0.48) in the positive CC-based graphs. This lack of sparsity in positive CC-based graph networks could be a result of the presence of the common noise process and global signal across the brain. Positive CCs were significantly affected by low-frequency physiological noise, which supports the notion that residual or unmodeled physiological noise (Bianciardi et al., 2009) creates a bias in internode CCs. RsFC connections characterized by negative CCs were not significantly affected by vasomotion arising from respiration and cardiac fluctuations. The sparse nature of negative CC graphs (presumably uncontaminated by common noise processes and global signal) indicates that anticorrelations in rsFC networks can be more specific to neuronal functional connectivity than correlations.
Methodological considerations
One of the drawbacks in this preliminary study is the loss of specificity engendered by replacement of hubness scores at a given voxel location by the maximum hubness value within 5 mm of that voxel. Thus, the location of the HACs that contributed to group hubness maps could vary by 5 mm across the subjects. However, the results from Table 1 support the adoption of this approximation in the methods. It also must be noted that the quantification of hubness based on DC can potentially inflate the importance of HACs, which are anticorrelated with voxels in brain networks with large spatial extent (Power et al., 2013). However, since in this study the HACs have been formulated to primarily highlight anticorrelated networks, the bias in location and centrality of the HACs does not detract from the goals of the project.
Furthermore, due to the proof-of-concept nature of this study, very conservative subject-specific brain GM masks were employed, through which voxels at brain edges and high susceptibility contrast boundaries were removed, thereby minimizing spurious sources of signal change. On the other hand, this conservative masking resulted in the loss of sensitivity to some brain voxels. Finally, in this study, we did not include the effects of dynamic changes in functional connectivity networks that have been reported recently (Sakoglu et al., 2010). Accounting for these dynamic changes in network behavior may better resolve different mechanisms of anticorrelations. The methodological shortcomings of this proof-of-concept study can be easily overcome in a future large sample project designed to better understand anticorrelated networks.
Conclusion
In this preliminary study, we found dense HACs in sensorimotor and emotion networks apart from basal ganglia and the thalamus. Examination of the anticorrelation maps of these hubs indicated that negative CCs in rsFC network architecture highlight a number of regulatory interactions between brain networks and regions, including reciprocal modulations, suppression, inhibition, and neurofeedback. Future research is needed to better understand and isolate anticorrelated networks in brain rsFC architecture.
Footnotes
Acknowledgments
This study was supported by the Department of Radiology and Imaging Sciences, Emory University and the Atlanta VA Medical Center.
Author Disclosure Statement
No competing financial interests exist.
