Abstract
Analyses of intrinsic fMRI BOLD signal fluctuations reliably reveal correlated and anticorrelated functional networks in the brain. Because the BOLD signal is an indirect measure of neuronal activity and anticorrelations can be introduced by preprocessing steps, such as global signal regression, the neurophysiological significance of correlated and anticorrelated BOLD fluctuations is a source of debate. Here, we address this question by examining the correspondence between the spatial organization of correlated BOLD fluctuations and correlated fluctuations in electrophysiological high γ power signals recorded directly from the cortical surface of 5 patients. We demonstrate that both positive and negative BOLD correlations have neurophysiological correlates reflected in fluctuations of spontaneous neuronal activity. Although applying global signal regression to BOLD signals results in some BOLD anticorrelations that are not apparent in the ECoG data, it enhances the neuronal-hemodynamic correspondence overall. Together, these findings provide support for the neurophysiological fidelity of BOLD correlations and anticorrelations.
Introduction
Examination of correlated spontaneous fluctuations in the BOLD signal represents an increasingly popular approach to the study of brain function and pathophysiology. Functional networks revealed by BOLD intrinsic functional connectivity (BOLD-iFC) are strikingly similar to those elicited with task-based paradigms (Smith et al., 2009; Laird et al., 2011; Mennes et al., 2012). However, our understanding of the neurophysiological significance of BOLD-iFC is complicated by non-neuronal contributions to the BOLD signal (e.g., cardiac and respiratory signals) (Birn et al., 2008; Chang and Glover, 2009b), motion artifacts (Power et al., 2012; Van Dijk et al., 2012), and image processing (Murphy et al., 2009; Weissenbacher et al., 2009). Because electrophysiological recordings of neuronal activity are less susceptible to these factors and are more directly interpretable, investigations of the electrophysiologic correlates of BOLD-iFC can help clarify its physiological bases (Scholvinck et al., 2010; Leopold and Maier, 2011).
Invasive electrode sampling for seizure localization provides an opportunity to directly record from the human cortex with electrocorticography. Local neuronal activity may be estimated from the electrocorticogram by calculating high (50–150 Hz) γ power (HGP) (Manning et al., 2009; Miller, 2010), which is the best known electrophysiological correlate of the evoked BOLD response (Logothetis et al., 2001; Mukamel et al., 2005; Conner et al., 2011). At rest, fluctuations in the HGP signal are correlated within auditory (Nir et al., 2008) and sensorimotor cortex (He et al., 2008). The correspondence between BOLD-iFC and HGP-defined functional networks demonstrated by these studies supports a neurophysiological basis for BOLD-iFC. However, a systematic exploration of the correlation structure of BOLD and HGP is required to determine the extent to which such neuronal-hemodynamic correspondence exists throughout cortex.
These investigations may prove particularly informative with regard to anticorrelated BOLD signal fluctuations, the electrophysiological correlates of which are unknown. Anticorrelated BOLD-iFC is a reproducible phenomenon (Shehzad et al., 2009) commonly observed between networks supporting apparently competing processes (Greicius et al., 2003; Fox et al., 2005), suggesting “functional segregation” (Fair et al., 2007; Kelly et al., 2008; Gee et al., 2011) possibly generated by (direct or indirect) inhibitory interactions (Greicius et al., 2003; Fox et al., 2012). However, interpretations of BOLD anticorrelations have been challenged because global signal regression (GSR), which removes a “global” BOLD signal to “unmask” physiologically significant relationships (Fox et al., 2009), also shifts the distribution of correlations leftward, possibly introducing spurious anticorrelations (Murphy et al., 2009; Weissenbacher et al., 2009).
Here, we hypothesized that the strength of correlated and anticorrelated BOLD signals would predict the strength of correlated and anticorrelated neuronal fluctuations measured directly from the cortical surface. We compared the correlation structure of BOLD and HGP fluctuations in five human subjects. After observing strong intraindividual correspondence throughout the cortex, we demonstrate neuronal anticorrelations in the human brain, which correspond spatially to regions of anticorrelated BOLD fluctuations. Finally, we demonstrate that GSR of BOLD data enhances the neuronal-hemodynamic correspondence and improves signal detection of underlying neuronal fluctuations. Together, these findings provide support for the neurophysiological fidelity of BOLD correlations and anticorrelations.
Materials and Methods
Patient selection.
Five patients (3 female, 33.8 ± 16.2 years of age; range, 22–60 years) with refractory epilepsy at Long Island Jewish Medical Center and North Shore University Hospital participated. All patients provided fully informed consent according to National Institutes of Health guidelines, as monitored by the local Institutional Review Board, in accordance with the ethical standards of the Declaration of Helsinki. Patient characteristics can be found in Table 1. To approximate normal brain physiology in this clinical population, subjects were selected out of 10 consecutive patients chosen for electrode implantation for seizure monitoring who had the following: (1) no prior resection, (2) brain MRI without extensive cortical dysplasias, tumors, or encephalomalacia, and (3) interictal epileptiform discharges fewer than one per minute throughout the duration of the recording session, which were removed before analysis. Electrocorticography (ECoG) was completed over the course of clinical monitoring for spontaneous seizures for localization of seizure foci. The decision to implant, the electrode targets, and the duration of implantation was made entirely on clinical grounds without reference to this investigation.
Clinical characteristics
Electrode implantation and recording.
Patients were implanted with intracranial subdural grids, strips, and/or depth electrodes (Integra Lifesciences) for 5–10 d. They were monitored in a specialized hospital setting until sufficient data were collected to identify the seizure focus. Continuous intracranial video-EEG monitoring was performed with standard recording systems (XLTEK EMU 128, LTM System) with a sampling rate of either 1000 or 500 Hz. A vertex screw into the bone was used as the reference for all electrodes. Electrophysiological data analyses were performed using custom MATLAB scripts (MathWorks).
Electrode registration.
To map electrophysiological findings to anatomical locations on the cortex, subdural electrodes were identified on the preoperative MRI by first registering the locations of the electrodes on the postimplantation CT to the equivalent location in the postimplantation structural MRI (Keller et al., 2011). Preimplantation and postimplantation MRIs were both skull-stripped using the BET algorithm from the FSL software library (www.fmrib.ox.ac.uk/fsl/) followed by coregistration to account for possible brain shift caused by electrode implantation and surgery (Mehta and Klein, 2010). Electrodes were identified in the CT using BioImagesuite (Duncan et al., 2004) and subsequently mapped to the closest point on the reconstructed pial surface of the preimplantation MRI in MATLAB using custom scripts (Dykstra et al., 2012). The reconstructed pial surface was computed using Freesurfer (Dale et al., 1999). Intraoperative photographs were used to corroborate this registration method based on the identification of major anatomical features.
Imaging.
Patients were scanned on a General Electric Signa HDx 3T scanner. Resting state fMRI data were acquired using one of two EPI gradient echo sequences (values in parenthesis represent parameters for patients 1, 4, and 5): FOV = 220 (240 mm), voxel size 4 × 3.5 × 3.5 (3.4 × 3.4 × 3.4 mm), matrix 64 × 64, flip angle = 50 (77), TR = 2000 ms, TE = 30, acquisition plane = axial, 150 contiguous volumes, 40 slices. Participants were instructed to rest with their eyes closed for 5 min. An anatomical T1-weighted image was acquired using one of two spoiled gradient-recalled acquisition in a steady state sequences: FOV = 256 mm (240 mm), voxel size 1×1×1 (1.2×0.9×0.9), matrix 256×256, flip angle = 8, TR = 7.8 ms (6.5 ms), TE = 3 (2.8), TI = 650 ms (600 ms), acquisition plane = axial (sagittal), 180 (170) slices.
BOLD intrinsic functional connectivity.
Resting state BOLD preprocessing was performed using AFNI (Cox, 1996) and FSL (Smith et al., 2004), and included slice timing correction for interleaved slice acquisition, volumetric realignment, despiking, spatial smoothing (6 mm full-width at half-maximum Gaussian blur), bandpass filtering (0.009–0.1 Hz), and linear and quadratic detrending. GSR was performed as part of a preprocessing step in which we regressed each patient's preprocessed time series on nine nuisance covariates (six head motion parameters, signals derived from the ventricles, white matter, and the global signal). To ensure that partial voluming of gray matter did not occur, we required 40% tissue type probability for the CSF mask and 66% tissue type probability for the white matter mask (Kelly et al., 2009). To examine the effects of GSR on the correspondence between BOLD and HGP correlations, resting BOLD data were also analyzed without GSR (i.e., each patient's preprocessed time series was regressed on eight nuisance covariates, comprising six head motion parameters, and signals derived from the ventricles and white matter). The resultant 4-D residuals volumes were registered to the patient's anatomical image using a linear transformation with 6 degrees of freedom.
For each patient and each recording site, spherical seed ROI (6 mm radius) were constructed centered in the gray matter located underneath each electrode position (“electrode-linked ROIs”). The mean time course for the seed was computed by averaging across all voxels within the ROI. A whole-brain map of BOLD-iFC was created by computing the correlation between the seed time course and that of every other voxel in the brain, and applying Fisher's r-to-z transformation to the resultant correlation coefficients. For visualization purposes, the resultant BOLD-iFC map was thresholded at |z| ≤ 0.3, projected from a 3D volume to the individual's pial surface, and plotted using MATLAB.
Electrocorticography.
ECoG was acquired for 3–6 min while subjects were asked to rest quietly. Two patients rested quietly with their eyes closed, whereas three rested with their eyes open. Interictal discharge free periods (276.1 ± 71.2 s, mean ± SD) were selected for analysis. Resting state sessions were conducted at least 2 h before or after an ictal event. To compare BOLD signals and electrophysiological measures on similar time scales, we computed band-limited power (BLP), a measure of the power fluctuations in specific frequency bands (Logothetis et al., 2001; Nir et al., 2008; Conner et al., 2011; Leopold and Maier, 2011; Ossandon et al., 2011). Channels with significant amounts of noise (SD > 250 uV) as well as electrode sites corresponding to the seizure onset zone were excluded (5.6 ± 3.2% of all channels, mean ± SD). The remaining channels were notch filtered (60 Hz) to remove power line noise and rereferenced by subtracting the common average to remove non-neuronal activity (Kanwisher et al., 1997). Data were bandpass filtered (fourth order zero-phase shift Butterworth filter) into equal 10 Hz frequency bands (e.g., 50–60, 60–70…140–150 Hz). For each 10 Hz band, the BLP signal was computed by full-wave rectification with the Hilbert transform to obtain the signal envelope (Crone et al., 2001; Nir et al., 2008; Conner et al., 2011; Leopold and Maier, 2011). BLP was then averaged across the 10 Hz bands to obtain the HGP. The band-specific normalization corrects for the 1/f drop-off of the power spectrum (Conner et al., 2011; Ossandon et al., 2011). Next, HGP signals were separated (fourth order zero-phase shift Butterworth filter) into slow (0.1–1 Hz) and fast (1–10 Hz) HGP fluctuations (Nir et al., 2008; Honey et al., 2012). HGP data <0.1 Hz were not assessed because of the small number of cycles recorded during the resting state experiment. Finally, the correlation between HGP signals (HGP-iFC) at each electrode was computed for each component of the HGP (slow, fast), and the resulting coefficients were transformed using Fisher's r-to-z transformation.
As the average reference reflects a type of normalization and may bias the findings presented here, we reanalyzed the ECoG by computing (1) a ground reference (no average reference) and (2) a local reference (Laplacian) wherein the average data from the electrodes immediately surrounding (<15 mm) the electrode of interest are subtracted before further analysis. In general, connectivity matrices of the HGP correlations computed with each method were very similar (rorig,avg = 0.91, rorig,laplacian = 0.90, ravg,laplacian = 0.94 for patient 4). Positive HGP correlations were stronger and more broadly distributed with the original ground reference compared with the average or Laplacian reference (data not shown). Negative HGP correlations, on the other hand, were located in similar locations and did not diminish in strength considerably.
GSR.
To approximate GSR as applied to BOLD-iFC, we computed the mean HGP from all nonartifactual electrodes (the “global” HGP signal) and regressed the signal at each electrode on the global HGP signal using linear regression. To assess the significance of the effect of GSR on the spatial correspondence of BOLD and HGP correlations, for each subject we first calculated the p value associated with the change in correspondence between BOLD and HGP correlations. Then, we performed a group analysis using Fisher's combined probability test, which combines independent tests from each subject that evaluates the same hypothesis (Fisher, 1925; Brown, 1975).
Identifying relationships between BOLD and HGP correlations.
To determine the association between the correlation structures of BOLD and HGP, HGP-iFCs were binned by BOLD-iFC strength, so that each bin contains HGP-iFC values whose complementary BOLD-iFC values are similar in strength (BOLD-iFC steps of r = 0.1 so bins are −0.2 < r <−0.1; −0.1 < r < 0; 0 < r < +0.1, etc). A scatterplot representation of these data before binning can be found in Figure 4 (data for patient 5). After HGP-iFC binning, an ANOVA was performed to test for the presence of a significant effect of BOLD-iFC strength on HGP-iFC. An ANOVA requires that data in each bin be normally distributed. Therefore, we assessed the normality of the HGP-iFC distribution in each bin by performing a Kolmogorov–Smirnov test. Any bins that were not normally distributed were removed from analysis. Next, a one-way ANOVA was performed, testing for a main effect of BOLD-iFC strength on HGP-iFC strength across a) all BOLD-iFC and b) negative BOLD-iFC only (see Fig. 3E). We examined the relationship between HGP-iFC and BOLD-iFC identified using each correlation technique (with and without GSR) between each electrode-linked ROI.
Assessing the statistical significance of negative HGP correlations.
As electrophysiological signals exhibit dependence on one another across electrodes, statistical significance of HGP correlations and anticorrelations was evaluated by a time-series block bootstrapping method (Hall and Horowitx, 1996; Politis, 2003). Briefly, for each electrode pair, one of the two time series of interest was broken into equal 10 s blocks and randomly shuffled in time. Block length was chosen to be large enough to maintain the temporal structure of the time series but small enough to permit several iterations of random shuffling of the dataset. The correlation coefficient between each randomly shuffled time series and nonshuffled time series was computed and repeated five times for each electrode pair, the average of which constituted one value in the null distribution. The location of each experimental correlation within the null distribution was assigned and a statistical cutoff for significance was determined using a false detection rate (FDR) to correct for multiple comparisons (q < 0.05) across electrode pairs (Benjamini and Hochberg, 1995; Benjamini et al., 2006).
Receiver operating curves.
To examine the spatial specificity of correlated and anticorrelated HGP-iFC to regions of correlated and anticorrelated BOLD-iFC, we computed receiver operating characteristic (ROC) curves for HGP-iFC (without GSR) and BOLD-iFC (with and without GSR). For each patient, a threshold for HGP-iFC was set as the 25th percentile of all negative HGP-iFCs. The number of anticorrelated HGP functional connections exceeding this threshold was then used to set a threshold for positive HGP-iFC, so that both sets of “target signals” (i.e., anticorrelated and correlated HGP-iFC) comprised an equal number of functional connections. Sensitivity and specificity with regard to correlated and anticorrelated BOLD-iFC were calculated using a several thresholds separated by 0.05 (e.g., r = 0, +0.05, +0.10; and r = 0, −0.05, −0.10).
Results
We investigated the electrophysiological correlates of correlated and anticorrelated BOLD fluctuations using direct intracranial recordings in five subjects (three female, age 22–60 years; see Table 1 for clinical characteristics) undergoing invasive intracranial evaluation with subdural electrodes. We identified networks of the following: (1) correlated and anticorrelated BOLD fluctuations (BOLD-iFC), and (2) correlated and anticorrelated HGP fluctuations (HGP-iFC) in preoperative resting state fMRI and postimplantation resting state ECoG data. An average of 100 electrodes were sampled in each patient (range, 51–122 electrodes), corresponding to 53,383 pairwise interactions across frontal, temporal, parietal, and occipital lobes.
Correlated and anticorrelated BOLD-iFC in patients with epilepsy
Epilepsy has been shown to alter resting state fMRI networks (Laufs et al., 2007; Broyd et al., 2009; Zhang et al., 2011). Therefore, we first verified that BOLD-iFC networks reported in healthy subjects were also observed in this patient population. Figure 1G depicts the voxelwise BOLD-iFC map of one patient, for a seed placed in dorsal prefrontal cortex, a prominent region within the default mode network (DMN). As expected, robust positive correlated BOLD-iFC was evident throughout the DMN, including in posterior cingulate cortex, lateral parietal cortex, and anterior middle temporal gyrus, whereas anticorrelated (negative) BOLD-iFC was observed in regions resembling intraparietal sulcus, frontal eye field, and supplemental motor area. DMN patterns were consistent across patients and visually exhibited similar spatial distributions with those obtained in healthy subjects (Fox et al., 2005).
Workflow. A, B, Coregistration of the preoperative MRI and postoperative CT scan with confirmation by intraoperative photograph allows electrode localization to anatomical features. B, Electrodes localized from A are superimposed on the MRI-reconstructed pial surface. C, D, Calculation of HGP. Raw voltage traces (B, green and black circles) were bandpass filtered (50–150 Hz) and enveloped to generate (D) slow (0.1–1 Hz) and fast (1–10 Hz) fluctuations of band-limited HGP. E, Correlation maps of HGP with seed at black circles are superimposed onto the pial surface. F, BOLD time series are extracted (B, from voxels under green and black circles) and preprocessed, including GSR. G, Voxelwise correlations of intrinsic BOLD fluctuations for a seed region in dorsal prefrontal cortex (black circle). Red and blue heat maps represent voxels exhibiting significant positive and negative BOLD correlations, respectively. There is significant positive BOLD-iFC at MPFC, posterior cingulate cortex (PCC), lateral parietal cortex (LP), and anterior middle temporal gyrus (aMTG) and negative correlations at intraparietal sulcus (IPS) and supplemental motor area (SMA).
Patterns of HGP-iFC exhibit strong correspondence with BOLD-iFC
To compare the patterns of correlated fluctuations in spontaneous BOLD activity and in the electrophysiological HGP signal, we produced spatial maps of BOLD-iFC and HGP-iFC for each electrode site. After coregistration and electrode localization (Fig. 1A,B), we derived BOLD-iFC by computing the correlation between the mean BOLD signal within electrode-linked ROIs (i.e., gray matter voxels falling within a spherical seed ROI underlying each electrode; Fig. 1; for details, see Materials and Methods and Keller et al., 2011), for all possible pairings of electrodes. Similarly, HGP-iFC was identified by converting voltage traces at each electrode to HGP and computing the pairwise correlations between HGP fluctuations at each electrode (Fig. 1C–E).
First, we investigated the spectral range of which the strength of BOLD-iFC predicts the strength of HGP-iFC recorded from the same subject. Across all subjects and networks sampled in the current report, the BOLD correlation structure exhibited significantly higher spatial correspondence with slow (0.1–1 Hz) HGP fluctuations than faster (1–10 Hz) HGP fluctuations (r = 0.40 compared with r = 0.28, p < 0.05, paired t test). As a result, the remainder of analyses is performed with slow (0.1–1 Hz) HGP signals and referred to as “HGP.”
To what degree are BOLD networks reflected in correlated fluctuations of neuronal electrical activity beyond primary sensory and motor cortex? As shown in Figure 2, correspondence between modalities was observed regardless of network sampled. Across patients and networks, the spatial correspondence between BOLD-iFC and HGP-iFC ranged from r = −0.05 to r = 0.77 (mean ± SD, 0.40 ± 0.12; Fig. 2). Electrodes with low correspondence may in part reflect measurement error introduced by coregistration, postimplantation brain shift, ROI size, as well as the nonsimultaneous collection of BOLD-iFCs and HGP-iFCs. As distance between sites likely modulates the strength of both neuronal and hemodynamic interactions between those sites (Salvador et al., 2005; Honey et al., 2009; Gee et al., 2011), a linear regression analysis was applied to determine the effect of distance on the correspondence between BOLD-iFC and HGP-iFC. At an exemplar region in one patient, the strength of the neuronal-hemodynamic correspondence was not diminished when interelectrode distance was included as a covariate (rpre = 0.58 to rpost = 0.69; data not shown). Across all electrodes and subjects, the mean correlation following linear regression of distance was reduced from rpre = 0.40 to rpost = 0.32.
Strong BOLD-HGP coupling across networks. A, Correspondence between HGP-iFC and BOLD-iFC associated with a single electrode site in patient 5 (seed region indicated by arrow). B, The color of each electrode represents the strength of the spatial correlation between BOLD-iFC and HGP-iFC at each electrode when the electrode of interest is used as a seed region.
To illustrate the correspondence across both positive and negative BOLD and HGP correlation, we identified exemplar networks in each of the five patients (Fig. 3A–D). As expected, performing GSR before the correlation analysis increased the number of regions exhibiting negative correlations and decreased the number of regions exhibiting positive correlations for both HGP and BOLD (Murphy et al., 2009). For each seed region, we observed the expected spatial overlap between regions exhibiting strong positive BOLD-iFC and those exhibiting positive HGP-iFC, both before and after GSR (Fig. 3A–D, black arrows). Regions exhibiting anticorrelated BOLD-iFC after GSR generally underlay electrodes exhibiting weakly positive or negative HGP-iFC (Fig. 3C, white arrows); in contrast, after GSR of HGP, the same electrodes exhibited strong negative HGP correlation (Fig. 3D, white arrows). Within each patient and across all cortical regions sampled (range 1–15 cm electrode separation), the strength of BOLD correlation significantly predicted the strength of HGP-iFC before GSR (Fig. 3E; one-way ANOVA, as an example for patient 1, Fig. 3E, top left: F(2,13) = 21.2, p < 0.001, Kolmogorov–Smirnov test of bin distribution, p < 0.05; statistical tests for patients 2–5 can be found in the lower left panels of Fig. 3E). GSR applied to BOLD signals did not change these relationships considerably (Fig. 3E, middle column). However, GSR applied to both HGP and BOLD signals increased the likelihood that regions showing anticorrelated BOLD-iFC would correspond to anticorrelated HGP-iFC, and weak BOLD-iFC would correspond to weak HGP-iFC (Fig. 3E, right column, one-way ANOVA after exclusion of positive BOLD-iFC, F(2,4) = 6.6, p < 0.001 for patient 1; patients 2–5 can be found in the lower right panels of Fig. 3E).
Spatial correspondence between HGP-iFC and BOLD-iFC. A, B, Voxelwise BOLD-iFC maps (A) before and (B) after GSR. Red and blue heat maps represent voxels with significant positive and negative BOLD-iFC, respectively, compared with the seed region (black circle). C, D, HGP-iFC map for five subjects using a single electrode seed (C) before and (D) after GSR. Black arrows represent regions of strong positive (A, B) BOLD-iFC and (C, D) HGP-iFC, whereas white arrows represent regions of (B) negative BOLD-iFC, (C) weak HGP-iFC, and (D) negative HGP-iFC. E, Whole brain comparison of all seeds with respect to BOLD-iFC and HGP-iFC before and after GSR of both modalities. Significance values represent one-way ANOVA test for the effect of BOLD-iFC on HGP-iFC for all data (upper values) and only negative BOLD-iFC (lower values). Each graph represents data from the patient corresponding to that row. Error bars indicate SE.
Effects of GSR on the neuronal-hemodynamic coupling
It has been suggested that GSR of resting state fMRI data is required to remove a confounding global signal, thus “unmasking” neurophysiologically valid relationships among large-scale networks (Fox et al., 2009). Figure 3 suggests that GSR may increase the neuronal-hemodynamic correspondence. To quantify this effect, we performed GSR before calculating both HGP-iFC and BOLD-iFC and evaluated the effect of GSR on the correspondence between modalities. Figure 4 illustrates the increase in the spatial correspondence between BOLD-iFC and HGP-iFC after GSR applied to both modalities in one patient (Fig. 4A,B; rpreGSR = 0.39, rpostGSR = 0.50, n = 2601 interactions, zdifference in correlation coefficients = −6.7, p < 0.001) and across patients (Fig. 4C). Relative to the neuronal-hemodynamic correspondence obtained in the absence of GSR for either modality, performing GSR on the BOLD data before correlation analysis increased the correspondence (rpreGSR = 0.35, rpostGSR = 0.40 for all patients; Fig. 4C). This increase was statistically significant in 4/5 patients, an effect that is significant at p = 0.00001, according to a two-tailed Fisher's combined probability test. Performing GSR on both modalities before correlation analysis further increased the correspondence (r = 0.41; Fig. 4C). This additional increase was statistically significant in 3 of 5 patients (p = 0.0001, two-tailed Fisher's test).
Effect of GSR on neuronal-hemodynamic correspondence. A, B, Whole brain spatial correspondence of HGP and BOLD-iFC (n = 2601 interactions) (A) before and (B) after GSR in one patient. C, The change in spatial correspondence for all interactions in all patients after performing GSR on BOLD-iFC and HGP-iFC. Each patient is represented by one continuous line. **p < 0.0001, Fisher's combined probability test.
Anticorrelated HGP fluctuations exist before signal regression
Next, we investigated the evidence for neuronal anticorrelations in the human brain. A time-series block bootstrapping technique was performed, producing a null distribution against which experimental HGP correlations were compared and from which correlation thresholds were derived (see Materials and Methods). For each subject, statistically significant anticorrelated HGP fluctuations (0.1–1 Hz; no regression) were observed between an average of 2.7% (range, 0.57–4.95%) of all possible electrode pairs (Fig. 5; mean 141 ROI pairs, range 56–323 ROI pairs per patient), whereas significant positive HGP correlations were observed between an average of 35.6% of electrode pairs (2070 ROI pairs). As expected, GSR applied to HGP increased the number of statistically significant anticorrelated HGP fluctuations to 30.6% (range, 19.09–39.8%) of all possible electrode pairs (Fig. 5; mean 1635 ROI pairs, range 418–2579 ROIs per patient) of all correlations, whereas significant positive HGP correlations were observed between an average of 24.1% of electrode pairs (1057 ROI pairs).
Significant HGP anticorrelations are present without GSR. Histograms of the distribution of shuffled and experimental data for HGP correlations before and after GSR are shown. The transparent layer outlines the distribution of the shuffled null dataset. Dotted red line indicates the significance threshold for anticorrelated HGP determined from the corresponding shuffled dataset. Values represent the percentage and absolute number of significant anticorrelated HGP regions out of all possible ROI interactions. *p < 0.05 after FDR correction.
Detection of positive HGP correlations using positive BOLD correlations
Given the observation of a positive correlation in BOLD-iFC, what is the probability that there exists an underlying neuronal correlation? How does GSR affect this probability? First, we provide an example of the spatial specificity of neuronal and BOLD interactions with a seed placed in medial prefrontal cortex (MPFC) in one patient. Before GSR, highly positive BOLD-iFC mapped to regions of highly positive HGP-iFC (Fig. 6A, black circle). Across all networks and patients and before GSR of BOLD-iFC, 28.2% of electrode pairs (162 neuronal correlations per patient) exhibiting positive BOLD-iFC also exhibited positive HGP-iFC, compared with 9.2% (66 neuronal correlations) of pairs exhibiting nonsignificant BOLD-iFC (p < 0.01, paired t test), and 8.1% (11 neuronal correlations) of pairs exhibiting anticorrelated BOLD fluctuations (Fig. 6B, p < 0.01). After GSR of BOLD data, 41.8% (149 neuronal correlations) of electrode pairs exhibiting positive BOLD correlation also exhibited positive HGP correlation, compared with 11.2% (78 correlations) of pairs exhibiting nonsignificant BOLD correlation (p < 0.01, paired t test), and 6.4% (13 correlations) of pairs exhibiting anticorrelated BOLD (Fig. 6B, p < 0.01). In summary, after GSR of BOLD data, we observed a stronger spatial correspondence between regions exhibiting positive neuronal correlations and those exhibiting positively correlated BOLD fluctuations.
Detection of correlated and anticorrelated HGP fluctuations using resting state BOLD correlations. A, Spatial correspondence of BOLD-iFC and HGP-iFC before GSR with seed in MPFC (black circle). There is strong positive BOLD- and HGP-iFC (black circles) and negative BOLD- and HGP-iFC (white boxes). B, C, ROC curves for detection of HGP-iFC using (B) correlated and (C) anticorrelated BOLD-iFC. For each subject, sensitivity and specificity measures were calculated before (left) and after (right) global regression of BOLD signals. Each patient is represented by a colored plot. BOLD-iFC thresholds are indicated by the black, white, and gray circles for each patient. The dotted line indicates the chance discrimination level. Subpanels, Percentage spatial overlap of BOLD-iFC and HGP-iFC with respect to the categories of positive, zero, and negative correlation. Numbers in bars represent the average number of (B) correlated and (C) anticorrelated HGP-iFC electrode pairs corresponding to BOLD-iFC voxel pairs of each category.
Next, we investigated the discriminability of positive BOLD-iFC with regard to positive HGP correlations by generating ROC curves. ROC analysis offers quantitative performance results from a binary classifier system as a discrimination threshold is varied. To define “true correlations,” we isolated the most positive HGP-iFC values (see Materials and Methods). We then constructed ROC curves assessing how well the true correlations could be detected by thresholding the BOLD-iFC values at various thresholds. For each threshold, specificity was defined as the proportion of HGP-iFC that are not positive and were identified as such with BOLD-iFC, whereas sensitivity was defined as the proportion of positive HGP-iFC, which are correctly identified with positive BOLD-iFC. Anticorrelated HGP signals were isolated and defined in a similar manner and are described later.
Both before and after GSR of BOLD data, positive BOLD-iFC exhibited higher than random discriminability (Fig. 6B, dotted line). GSR appeared to reduce the intersubject variability in discriminability (Fig. 6B, compare left and right panels). Additionally, performing GSR on BOLD signals increased the overall discriminability (area under the curve). At a weak threshold (r > 0.1) for correlated BOLD-iFC, mean sensitivity and specificity were 90.0% and 27.7% before GSR of BOLD-iFC and 83.7% and 59.1% after GSR, respectively (Fig. 6B, gray circles). At a more stringent threshold (BOLD-iFC, r > 0.3), mean sensitivity/specificity was 64.2%/72.8% and 53.4%/91.3% before and after GSR, respectively (Fig. 6B, white circles). This analysis demonstrates that, overall, regions of positively correlated BOLD-iFC reliably discriminate positive neuronal correlations, and this discriminability increases with the application of GSR to BOLD signals.
Detection of negative HGP correlations using negative BOLD correlations
Given the observation of an anticorrelation in BOLD-iFC, what is the probability that there exists an underlying neuronal anticorrelation? For the seed placed in MPFC in one patient, regions of anticorrelated BOLD-iFC corresponded relatively well to regions of anticorrelated HGP-iFC (Fig. 6A, white boxes). Across all networks and patients and before GSR of BOLD data, 25.1% of electrode pairs (87 neuronal anticorrelations) exhibiting anticorrelated BOLD-iFC also exhibited anticorrelated HGP-iFC, compared with 10.1% (120 anticorrelations) of pairs exhibiting nonsignificant BOLD-iFC (p < 0.05, paired t test), and 11.4% (104 anticorrelations) of pairs exhibiting positive BOLD-iFC (Fig. 6B, p < 0.05). After GSR of BOLD-iFC, 21.9% of electrode pairs (217 neuronal anticorrelations) exhibiting anticorrelated BOLD-iFC also exhibited anticorrelated HGP-iFC, compared with 12.8% (101 anticorrelations) of pairs exhibiting nonsignificant BOLD-iFC (p < 0.01, paired t test), and 4.2% (43 anticorrelations) of pairs exhibiting positive BOLD-iFC (Fig. 6B, p < 0.01). In summary, after GSR of BOLD data, we observed the following: (1) an increase in the number of neuronal anticorrelations in regions of anticorrelated BOLD-iFC accompanied by (2) improved spatial specificity of neuronal anticorrelations to regions of anticorrelated BOLD-iFC, relative to regions of positive BOLD-iFC.
With regard to the ability of anticorrelated BOLD-iFC to discriminate neuronal anticorrelations, above-chance discriminability was observed for all 5 patients (Fig. 6C, dotted lines). As was the case for positive correlations (i.e., Fig. 6B), GSR applied to the BOLD signals increased the area under the ROC curve (Fig. 6C, compare left and right panels). Across subjects, sensitivity reached 46.4% before GSR and 76.8% after GSR (at r < 0 BOLD-iFC threshold). Without GSR, sensitivity and specificity measurements were to 28.8% and 86.4%, respectively, at a stronger BOLD-iFC threshold (r < −0.1) but were 60.0% and 65.6%, respectively, after GSR (Fig. 6C, gray circles). Sensitivity and specificity measurements at a higher BOLD-iFC threshold (r < −0.2) after GSR were 41.5% and 78.5%. These data demonstrate that, overall, neuronal anticorrelations are reliably indicated by anticorrelated BOLD-iFC after GSR.
Discussion
Here, we show the following: (1) correspondence between HGP-iFC and BOLD-iFC is robust across a variety of functional networks, extending previous studies; (2) a subset of networks exhibit significant neuronal anticorrelations, and these tend to correspond to BOLD anticorrelations; (3) neuronal correlations and anticorrelations exhibit strong spatial overlap with correlated and anticorrelated BOLD fluctuations, although the overlap is stronger for positive correlations; and (4) performing GSR on BOLD data enhanced this correspondence but also introduced anticorrelations that lacked obvious electrophysiological counterparts.
From BOLD to neural activity
HGP activity has been shown to correlate with neuronal firing rates (Manning et al., 2009), and changes in the evoked BOLD signal are best explained by changes in HGP (Logothetis et al., 2001; Mukamel et al., 2005; Conner et al., 2011). Therefore, HGP in awake humans appears to provide the closest in vivo measure of large-scale neuronal dynamics (Miller, 2010). If HGP represents the electrophysiological correlate of the BOLD signal, spontaneous fluctuations in HGP should similarly reflect spontaneous fluctuations in the BOLD signal. We observed robust spatial overlap between functional networks defined by correlated HGP fluctuations (i.e., HGP-iFC) and correlated BOLD fluctuations (i.e., BOLD-iFC). These findings constitute an important extension of previous reports of the correspondence between HGP-iFC and BOLD-iFC within visual cortex in primates (Scholvinck et al., 2010), and within auditory (Nir et al., 2008) and sensorimotor (He et al., 2008) cortex in humans, by demonstrating intraindividual correspondence in multiple networks supporting diverse cognitive functions (e.g., the default network) as well as networks supporting sensory and motor function. Our findings thus suggest that the neuronal-hemodynamic correspondence is not limited to specialized networks but constitutes a physiological property of functional networks throughout the brain.
Although a small subset of networks exhibited weak HGP-iFC/BOLD-iFC correspondence, no consistent pattern that accounted for the locations of such networks was apparent. Potential explanations for networks exhibiting weak correspondence include electrode coregistration errors, brain shift after implantation, seed ROI size, and the fact that the recording sessions were not simultaneous. In contrast to the minority of networks exhibiting weak neuronal-hemodynamic correspondence, the majority of networks exhibited relatively strong correspondence (rmean = 0.40, rmax = 0.77), suggesting that, on average, between 16% (rmean2 = 0.16) and 59% (rmax2 = 0.59) of the variance of BOLD-iFC can be explained by HGP-iFC recorded during a separate resting session. If simultaneous recordings of BOLD and intracranial electrophysiological data were possible in humans, we expect that the correspondence of these intrinsic signals would increase significantly. In previous studies, the percentage variance explained by neuronal signals has reached 10% during simultaneous monkey recordings (Scholvinck et al., 2010), 34% in the sensorimotor cortex in humans (He et al., 2008), and 50% under optimal conditions of simultaneous recordings of stimulus-driven responses (Logothetis, 2002). Although our analysis of neuronal activity is restricted to HGP, it is possible that band-limited power fluctuations in other frequencies (i.e., α, β) also show good correspondence with BOLD-iFC. We focused on HGP because it is thought to be most closely related to neuronal firing (Manning et al., 2009; Miller, 2010). However, given that HGP also exhibits modulation by the phase of lower-frequency bands (Lakatos et al., 2008), it is possible that HGP fluctuations may be mediated by lower-frequency phenomena. For example, the slow cortical potential (<.5 Hz) has been shown to correlate with BOLD fluctuations at least for local interactions in the sensorimotor cortex (He et al., 2008). Our preliminary data suggest that the slow cortical potential does partially reflect BOLD-iFC but that HGP fluctuations demonstrate stronger spatial correspondence with BOLD-iFC, especially when examining long-range interactions (data not shown). However, this question warrants further examination with a systematic analysis in follow-up studies.
Anticorrelated neuronal fluctuations
Here, we provide direct evidence of anticorrelated intrinsic neuronal activity in humans at rest. Although the strength of neuronal anticorrelations was lower than that of positive neuronal correlations, a subset did exceed a conservative FDR threshold for statistical significance. Given the local nature of neuronal activity in the high γ band (Crone et al., 2001; Miller et al., 2007), it is unlikely that electrical interference from neighboring cortical regions can account for these findings. Although anticorrelations have been observed in BOLD-iFC before GSR (Chang and Glover, 2009a; Chai et al., 2012; Liang et al., 2012), these findings represent the first report of anticorrelated neuronal activity in resting humans.
Time periods characterized by anticorrelated HGP fluctuations have been observed in electrophysiological recordings from task-on and task-off regions of cat cortex (Popa et al., 2009). In contrast, MEG studies examining band-limited power in humans have failed to observe anticorrelations in their recordings (de Pasquale et al., 2012). This divergence in findings may be reflective of differences in the frequency band examined (β vs high γ) or the recording technique (MEG vs EEG). We previously investigated the neurophysiological underpinnings of BOLD-iFC using corticocortical-evoked potentials (CCEPs) elicited by direct cortical stimulation in similar subjects. Although regions of positive BOLD-iFC corresponded to regions exhibiting a high degree of connectivity defined by CCEPs, no consistent CCEP correlate was seen for anticorrelated BOLD fluctuations (Keller et al., 2011). However, CCEPs probe the brain's response to exogenous input (its effective connectivity). In contrast, HGP-iFC recorded from patients in the resting state reflects a spontaneous functional connectivity and may be a closer analog of resting BOLD-iFC.
The spontaneous neuronal anticorrelations observed here likely reflect the electrophysiological correlate of anticorrelated BOLD signals. Electrodes exhibiting anticorrelated HGP-iFC were distributed across regions of correlated, noncorrelated, and anticorrelated BOLD-iFC but were significantly more likely to occur in regions of anticorrelated BOLD-iFC. The observation of neuronal anticorrelations outside regions of anticorrelated BOLD-iFC is perhaps not surprising as the neural-hemodynamic relationship is not a 1:1 relationship. In support of this, a recent study demonstrated that neuronal deactivation during goal-directed tasks, thought to be confined to the DMN when measured as BOLD signal decreases (Miller et al., 2009; Ossandon et al., 2011), was also observed outside this network (Ramot et al., 2012). Importantly, we found that the sensitivity of anticorrelated BOLD-iFC after GSR with regard to neuronal anticorrelations reached 76.8%. Although the discriminability (area under the ROC curve) and sensitivity of neuronal correlations in positive BOLD-iFC exceeded that of anticorrelated BOLD-iFC (up to 83.7%), both exhibited notable signal detection capabilities. Together, these findings support the neurophysiological veracity of anticorrelated BOLD-iFC.
GSR: revealing or obscuring?
The application of GSR to BOLD during the computation of iFC continues to be actively debated for several reasons. The primary motivation for GSR is the removal of non-neuronal noise generated by a variety of sources, such as instability of MRI measurements, heart rate, and respiration (Wise et al., 2004; Bianciardi et al., 2009), which can help unmask “true” interregional relationships (Fox et al., 2009). Although GSR will redistribute correlations to be zero-centered (Murphy et al., 2009), the spatially specific, reproducible anticorrelations it reveals appear to constitute neurophysiologically relevant relationships among regions and networks (Fox et al., 2009; Shehzad et al., 2009). The global BOLD signal is distributed heterogeneously throughout the cortex (Fox et al., 2009). Thus, GSR not only induces spurious negative correlations (Murphy et al., 2009; Weissenbacher et al., 2009) but may alter the strength and spatial distribution of positive correlations (Saad et al., 2012). Importantly, GSR may also remove a global neuronal signal arising from a widely distributed ascending input (Fox et al., 2009; Scholvinck et al., 2010), possibly reflecting a true shared covariation in firing rate. If so, this signal would be justifiably removed by GSR to reveal neuronal relationships otherwise masked by the dominant global signal (Fox et al., 2009).
The present results indicate that, although GSR can introduce some artifacts, it likely reveals more than it obscures. GSR significantly improved the correspondence between BOLD-iFC and HGP-iFC. Performing GSR before computing HGP-iFC increased the correspondence further, potentially reflecting the benefits of removing an obscuring global signal from both the neuronal and BOLD data. Furthermore, when evaluating the accuracy of BOLD-iFC for detection of HGP-iFC, GSR performed on BOLD-iFC decreased the variability in discriminability between subjects while increasing the sensitivity and specificity of BOLD-iFC with regard to both correlated and anticorrelated HGP-iFC.
It is important to underscore that the improved neuronal-hemodynamic correspondence of spontaneous fluctuations with GSR is not without potential costs or caveats. As expected, GSR shifted the distribution of correlations to be near zero-centered (Murphy et al., 2009), potentially leading to artifactual negative correlations, consistent with the predictions of critics of GSR. Although it is tempting to argue that the increased neuronal-hemodynamic correspondence of GSR justifies or vindicates the methodology, it is possible for an artifactual method to erroneously pull two distributions closer to one another. Moreover, although these results suggest that GSR enhances the neuronal-hemodynamic correspondence, they do not necessarily demonstrate that GSR is the most effective method of removing non-neuronal contributions to the BOLD signal. Partial correlation has been suggested as an alternative to GSR because it permits estimation of the relationship between two time series while controlling for signals from other sources of interest (Smith et al., 2011); importantly, the approach does not require an explicit modeling of a “global” signal. It is conceivable that partial correlation improves the neuronal-hemodynamic correspondence and detects underlying neuronal events as well or better than GSR. Unfortunately, partial correlation analysis could not be performed with our BOLD data as >100 regressors would be required. Longer scanning times, which increase the degrees of freedom available, may permit partial correlation analysis (Smith et al., 2012). Other alternatives to GSR, such as CompCor (Behzadi et al., 2007; Chai et al., 2012), which implements a component method to reduce noise in resting fMRI, allow for less problematic interpretations of BOLD-iFC. Understanding the effect of these alternatives on intermodal correspondence remains an important future endeavor.
Footnotes
This work was supported by the Page and Otto Marx Jr Foundation to A.D.M., the National Institute of Neurological Disorders and Stroke Grant F31NS080357-01 to C.J.K., the National Institute of Neurological Disorders and Stroke Grants F31NS080357-01 and T32-GM007288 to C.J.K., and the Epilepsy Foundation of America Grant EFA189045 to C.J.K. We thank Peter Kingsley for assistance with resting state fMRI sessions; Ido Davidesco, Meir Meshulam, and Rafael Malach for helpful discussions regarding BLP calculations; Miklos Argyelan, Gad Klein, and Andrew Dykstra for help with developing methodology to coregister electrode locations with noninvasive data; and the patients who participated in this study and the nursing and physician staff at North Shore–Long Island Jewish Medical Center.
The authors declare no competing financial interests.
- Correspondence should be addressed to Dr. Ashesh D. Mehta, Department of Neurosurgery, Hofstra North Shore LIJ School of Medicine and Feinstein Institute for Medical Research, 300 Community Drive, Manhasset, NY 11030. amehta{at}nshs.edu
