Abstract
Spontaneous activity observed with resting-state fMRI is used widely to uncover the brain's intrinsic functional networks in health and disease. Although many networks appear modular and specific, global and nonspecific fMRI fluctuations also exist and both pose a challenge and present an opportunity for characterizing and understanding brain networks. Here, we used a multimodal approach to investigate the neural correlates to the global fMRI signal in the resting state. Like fMRI, resting-state power fluctuations of broadband and arrhythmic, or scale-free, macaque electrocorticography and human magnetoencephalography activity were correlated globally. The power fluctuations of scale-free human electroencephalography (EEG) were coupled with the global component of simultaneously acquired resting-state fMRI, with the global hemodynamic change lagging the broadband spectral change of EEG by ∼5 s. The levels of global and nonspecific fluctuation and synchronization in scale-free population activity also varied across and depended on arousal states. Together, these results suggest that the neural origin of global resting-state fMRI activity is the broadband power fluctuation in scale-free population activity observable with macroscopic electrical or magnetic recordings. Moreover, the global fluctuation in neurophysiological and hemodynamic activity is likely modulated through diffuse neuromodulation pathways that govern arousal states and vigilance levels.
SIGNIFICANCE STATEMENT This study provides new insights into the neural origin of resting-state fMRI. Results demonstrate that the broadband power fluctuation of scale-free electrophysiology is globally synchronized and directly coupled with the global component of spontaneous fMRI signals, in contrast to modularly synchronized fluctuations in oscillatory neural activity. These findings lead to a new hypothesis that scale-free and oscillatory neural processes account for global and modular patterns of functional connectivity observed with resting-state fMRI, respectively.
Introduction
fMRI has been used increasingly to uncover large-scale neural networks with correlated spontaneous activity in the absence of any overt task (Biswal et al., 1995). These resting-state networks (RSNs) mostly appear modular (Van Dijk et al., 2010), arise from structural connections (Greicius et al., 2009), persist across behavioral states (Horovitz et al., 2008), and resemble task activation patterns (Smith et al., 2009), thereby being recognized to report on the brain's intrinsic functional organization (Fox and Raichle, 2007; Power et al., 2011; Yeo et al., 2011).
However, modular RSNs may be obscured by nonspecific global fluctuations in spontaneous fMRI activity (Fox et al., 2009). In practice, the global fMRI signal is often regressed out in preprocessing, but nevertheless creates confounding effects (Murphy et al., 2009; Saad et al., 2012) or alters clinical inferences (Yang et al., 2014). The origins of the global fMRI signal are elusive, with previous studies emphasizing the contributions from ongoing neuronal (Schölvinck et al., 2010) or non-neuronal (Birn et al., 2006; Shmueli et al., 2007) processes. How to characterize, distinguish, and interpret global versus modular RSN remains a lingering and fundamental challenge because fMRI measures the indirect hemodynamic signature of neural activity (Logothetis, 2008).
To address this issue, we used neurophysiological signals in multiple macroscopic scales, including electrocorticography (ECoG), magnetoencephalography (MEG), and electroencephalography (EEG). The common spectral, temporal, and spatial characteristics shared across such multilevel electrophysiological signals were identified and compared with those of resting-state fMRI. Our emphasis was on testing the hypothesis that global resting-state fMRI activity has a neural origin that manifests itself as the broadband power fluctuation of scale-free neurophysiology observable on the cortical or even head surface. Several previous findings led us to raise this hypothesis. Widespread fMRI activity was found to correlate with the powers of local field potential (LFP) consistently across a broad range of frequencies (>40 Hz; Schölvinck et al., 2010). Cortical electrical activity was predominantly arrhythmic and scale free in various arousal states (He et al., 2010). Scale-free neurophysiological signals were well characterized by power–law spectral distributions, for which changes in the broadband power of LFP, ECoG, and EEG, were correlated consistently with neuronal spiking activity during tasks or at rest (Manning et al., 2009; Miller et al., 2009; Whittingstall and Logothesis, 2009; Ray and Maunsell et al., 2011).
To test this hypothesis, we also used a recently developed computational tool (Wen and Liu, 2016) to separate scale-free and oscillatory components of neural signals so that the fluctuation and correlation of scale-free activity were assessed without any interference from brain rhythms and vice versa. Our results suggest that scale-free and oscillatory neural processes contribute differentially to resting-state fMRI: the former underlies the global signal, whereas the latter supports modular networks.
Materials and Methods
Data acquisition
EEG-fMRI.
EEG and fMRI data were acquired simultaneously from 19 healthy volunteers (age 28 ± 10, 9 female) in the eyes-closed resting state for 10 minutes. The EEG data were collected from 31 scalp channels (together with one unipolar ECG channel) with the standard montage and an MR-compatible recording system (BrainAmp MR Plus; BrainProducts). The EEG signals were referenced to the FCz electrode and sampled continuously at 5 kHz with a resolution of 0.5 μV/bit for a 16-bit amplifier and an analog bandwidth from 0.1 to 250 Hz. The EEG sampling clock was synchronized with an external reference signal obtained from the 10 MHz master clock of the MRI scanner. A slice-trigger signal that marked the onset time of every fMRI slice acquisition was recorded based on a TTL signal from the scanner. The fMRI acquisition was evenly spaced with equal delay between each slice acquisition.
The MRI and fMRI images were collected from a 3 tesla Signa MRI system (General Electric) equipped with a 16-channel receiver-only phase-array coil (NOVA Medical). Whole-brain fMRI data were acquired using a single-shot gradient recalled echo-planar imaging (EPI) sequence with 75° flip angle (FA), 30 ms echo time (TE), 1.5 s repetition time (TR), 30 axial slices with 4 mm thickness, 220 × 165 mm2 field of view (FOV) for 64 × 48 matrix size, and sensitivity encoding (SENSE) parallel imaging with an acceleration factor of two. The last imaging volume was acquired with slightly longer TE to measure the spatial variation of the B0 field. T1-weighted anatomical images covering the whole head were acquired using a 3D magnetization prepared rapid gradient echo (MPRAGE) sequence with 12° FA, 2.25 ms TE, 5 ms TR, 725 ms inversion time, 200 sagittal slices, and 1 mm3 isotropic resolution.
MEG.
MEG data were a part of the data from our previous study (Liu et al., 2010). Data from three human volunteers were used for the purposes of this study. Each subject was instructed to lie in a dark and magnetically shielded room. The experimental paradigm started with two cycles of alternating eyes-closed and eyes-open wakeful resting states and then the subjects fell asleep naturally. The total recording period including all three states was 40–45 min. The MEG data were recorded from 275 radial first-order gradiometer channels covering the entire head with a CTF MEG system. Signals were sampled at 600 Hz and background noise was removed online.
ECoG.
Macaque ECoG data were publicly available on the website of Project Tycho (http://neurotycho.org). These data were originally acquired from 128 sensors covering the entire lateral surface of the left macaque cortex with a sampling rate of 1 kHz using a Cerebus recording system (Blackrock Microsystems). Data were referenced to an electrode between the silicone sheet of the subdural electrodes and the dura mater with the platinum side facing the dura and the ground electrode in between the dura and skull with the platinum side facing the skull (for details, see Nagasaka et al., 2011). In this study, we downloaded and selected the data from two macaques (denoted as macaque C and G) in three behavioral states: eyes-closed wakefulness, eyes-open wakefulness, and sleep. For macaque C, there were ECoG data from six sessions in eyes-closed wakefulness, three sessions in eyes-open wakefulness, and three sessions in sleep. For macaque G, there were five sessions in eyes-closed wakefulness, two sessions in eyes-open wakefulness, and three sessions in sleep. Every session contained 5 or 10 min ECoG time series from 128 channels.
Data preprocessing
EEG.
The recorded EEG data contained gradient and cardiac ballistic artifacts from concurrent fMRI acquisition. Such artifacts were removed using software developed in-house (Liu et al., 2012a) and publicly available (https://purr.purdue.edu/publications/1936). After down-sampling the data to 1 kHz, a third-order polynomial function was used to model and fit the slow trend in EEG time series and then removed channel by channel. The detrended data were further low-pass filtered with a 250 Hz cutoff frequency.
fMRI.
The EPI images were corrected for geometric distortion using the B0 field map. The EPI time series were further corrected for slice timing differences and subtle head motion using FSL (FMRIB, Oxford, UK) and then were detrended by regressing out the third-order polynomial function that modeled the slow signal drift. A rigid-body transformation was estimated and used to align the EPI images to the T1-weighted MPRAGE images using align_epi_anat.py in AFNI (National Institute of Mental Health–National Institutes of Health). The T1-weighted images were transformed into the standard MNI space using a nonlinear registration tool (FNIRT) in FSL. Through this nonlinear transformation, the EPI data were normalized to the MNI space and resampled to 3 mm isotropic resolution. The normalized EPI data were spatially smoothed with a 3D Gaussian kernel with 4 mm full width at half maximum. Finally, every voxel time series was demeaned and standardized.
MEG.
Physiological noises of cardiac or respiratory origin were removed from MEG raw data using independent component analysis, as described in our previous study (Liu et al., 2010). Specifically, the components originating from cardiac or respiratory activity were identified in a semiautomatic manner. An autocorrelation function was computed based on the time course of each component. A component was identified as a cardiac artifact if a positive peak autocorrelation coefficient >0.3 was found at a time lag between 0.6 and 1.5 s or as a respiratory artifact if a positive peak correlation coefficient >0.3 was found at a time lag between 2.5 and 5 s. These particular intervals were chosen to cover the expected ranges of the cardiac and respiratory rate. The identified artifact components were inspected visually and confirmed and then excluded. The remaining components were reassembled to yield the denoised MEG data, and then the MEG signals were low-pass filtered with a 250 Hz cutoff frequency. As described previously and used elsewhere (Horovitz et al., 2008; Liu et al., 2010), the sleep period was determined using the inverse index of wakefulness (IIoW), defined as the ratio of the power in the delta and theta bands to that of the alpha band computed over 2 min intervals. This was because transition from wakefulness to light sleep was characterized by an attenuation in the alpha rhythm and an elevation in the theta and delta rhythms in both EEG and MEG (Simon et al., 2000; Olbrich et al., 2009). The sleep period was initially determined as a sustained period (>20 min) in which IIoW was twice as much as that during eyes-closed wakefulness and this was inspected visually and confirmed by the authors.
ECoG.
A third-order polynomial trend in ECoG time series was removed channel by channel. The detrended data were further low-pass filtered with a 250 Hz cutoff frequency.
Data analysis
Separating the power spectrograms of scale-free and oscillatory electrophysiology.
For ECoG, MEG, and EEG, we used the irregular-resampling auto-spectral analysis (IRASA) to separate the power spectrograms of underlying scale-free and oscillatory electrophysiology (Wen and Liu, 2016). A sliding time window was defined to be 3 s long and step by 1 s for ECoG and MEG and 0.125 s for EEG. For every time window, the preprocessed ECoG, MEG, or EEG time series was separated into two components reflecting scale-free and oscillatory neural dynamics. The scale-free component is arrhythmic and not confined to any specific time scale (Yamamoto et al., 1991); its power spectrum follows a power-law distribution, shown as a descending line on the log–log plot (Miller et al., 2009; He, 2014). In contrast, the oscillatory component contains rhythms indicative of neural dynamics with characteristic time scales (Buzsáki and Draguhn, 2004); it typically exhibits one or multiple spectral peaks at specific frequencies. Based on their distinct temporal and spectral characteristics, we were able to separate the scale-free and oscillatory components in the power spectrum of the neurophysiological signal.
Briefly, we irregularly resampled a neural signal by multiple pairs of non-integer factors. For each pair, one was chosen as a value between 1.1 and 1.9 with 0.05 increments and the other was chosen to be its reciprocal. The former was used to up-sample the signal, whereas the latter was used to down-sample the signal by the same factor. We then computed the geometric mean of the auto-power spectra of every pair of the resampled (i.e., up-sampled and down-sampled) signals. In the resulting spectrum, the power associated with the oscillatory component was redistributed away from its original (fundamental and harmonic) frequencies by a frequency offset that varied with the resampling factor, whereas the power solely attributed to the scale-free component remained the same power-law distribution independent of the resampling factor. By taking the median of the auto-power spectra of the variously resampled signals, we extracted the power spectrum of the scale-free component; the difference between the original power spectrum and the extracted broadband scale-free spectrum was an estimate of the power spectrum of the oscillatory component (see Fig. 1 for an example). Detailed procedures were described in our previous publication and tested with computer simulation and experimental data (Wen and Liu, 2016). Applying IRASA to electrophysiological signals within every sliding window, we obtained the temporally fluctuating scale-free and oscillatory power spectra; that is, time–frequency representations or spectrograms.
Broadband and narrow-band power fluctuations.
From the obtained spectrogram of the scale-free component of ECoG, MEG, or EEG, we extracted the temporal fluctuation of the broadband power for every recording channel. Specifically, for a given channel and a given time window, the power spectrum of the scale-free component was transformed to the double logarithmic scale. The scale-free spectrum did not exhibit as a strict line in log–log scale, but as two lines that joined at a “knee” frequency around 15 Hz. Therefore, we separated the frequency points into two ranges: one for the low frequencies (LFs) from 1 to 15 Hz and the other for the high frequencies (HFs) from 15 to 100 Hz. In each of these two ranges, the frequency points were resampled to be evenly distributed in terms of the log frequency; the log power was averaged across the resampled frequency points, yielding the average log power as a measure of the broadband power (BBP) in the corresponding frequency range. The BBP time series indicates the power fluctuation of frequency-independent scale-free electrophysiology. These procedures were also used in previous studies by us (Wen and Liu, 2016) and others (Miller et al., 2009).
From the separated spectrogram of the oscillatory component of ECoG, MEG, or EEG, we extracted the temporal fluctuation of the band-limited power for the alpha band (8–13 Hz). In this study, we only focused on the alpha-band-limited power (BLP) because a spectral peak in the alpha band was observed consistently in our datasets (Liu et al., 2012b). The BLP time series indicates the power fluctuation of frequency-specific oscillatory electrophysiology.
We also examined the relationship between the power fluctuations of scale-free electrophysiological signals in the LF and HF ranges. For this purpose, the cross-correlation between LF-BBP and HF-BBP was calculated channel by channel for ECoG and MEG. For every channel, the significance of the correlation was assessed by calculating the p-value with the degree of freedom (DOF) equal to the number of nonoverlapping time windows (in the time–frequency analysis) − 2 (DOF = 98 or 198 for the 5 min or 10 min ECoG or DOF = 398 for MEG with Bonferroni correction to account for the number of channels and p < 0.001). The percentage of significantly correlated channels was further calculated for each session and each arousal state. Note that, for the human MEG data, we only used the sleep period for calculating the cross-correlation such that the correlations were only attributed to activity fluctuations within a single arousal state instead of being dominated by activity differences among different states.
We also addressed whether the level of fluctuation in scale-free electrophysiology differed between different arousal states. For each ECoG or MEG channel, the fluctuation level of BBP was quantified as the ratio of the temporal SD to the temporal average. The difference in the fluctuation level between eyes-closed wakefulness and sleep was evaluated and tested for significance. Specifically, the BBP time course at each channel was first divided by its temporal average. For such “normalized” signals, a two-sample F test was used to test the statistical significance of the difference in variance between eyes-closed wakefulness and sleep (p < 0.01 with Bonferroni correction to account for the total number of channels). To sum across channels, we further calculated the percentage of the channels with significant between-state differences in the BBP fluctuation level.
Correlation in scale-free versus oscillatory ECoG and MEG.
After obtaining the temporal fluctuations of BBP (i.e., scale-free activity) and BLP (i.e., oscillatory activity) for every ECoG or MEG channel, we analyzed the spatial patterns of temporal correlations across channels. Such cross-correlations were analyzed for each pair of different channels and shown as a correlation matrix in which each row corresponds to a map of correlations to a seed channel. Example correlation maps with specific seed channels were displayed for illustration. The statistical significance of the cross-channel correlation was assessed with the DOF as the number of nonoverlapping time windows − 2 (DOF = 98 or 198 for ECoG, or DOF = 398 for MEG) and Bonferroni correction for multiple comparisons over all different pairs of channels with p < 0.001.
We further calculated the percentage of significantly correlated pairs of channels and assessed the inhomogeneity of pairwise correlations as the ratio of their SD to their average. For the ECoG data, these metrics were calculated separately for each arousal state (i.e., eyes-open and eyes-closed wakefulness and sleep). For the MEG data, they were only calculated for the sleep state because the period in either the eyes-open or eyes-closed state was too short (2 min) to yield reliable correlation measures. Note that a higher percentage of significantly correlated channels indicates a greater degree of global correlations. A lower level of inhomogeneity in such correlations suggests a more nonspecific (as opposed to modular) correlational structure. Altogether, these two metrics served to quantify the degree of nonspecific and global correlation observed with either scale-free or oscillatory electrophysiology.
Further, we investigated whether the degree of global and nonspecific correlation differed between the scale-free (i.e., BBP) and oscillatory (i.e., BLP) components of electrophysiology or across different arousal states. The percentage or inhomogeneity metrics based on BBP (i.e., scale-free) were compared against those based on BLP (i.e., oscillatory) across different sessions in the same state. Similarly, these metrics observed with the ECoG BBP were also compared between awake and sleep states. Based on the observation from the exploratory seed-based correlation analysis, we hypothesized that the degree of global and nonspecific correlation was higher in the scale-free component compared with the oscillatory component and during sleep compared with during wakefulness for the scale-free component. To test both hypotheses, one-tailed nonparametric Wilcoxon signed-rank test (or t test) was used to evaluate the significance of the difference in the percentage (or inhomogeneity) with p < 0.05. A nonparametric test was used because these variables were non-Gaussian or bounded.
In addition, we evaluated the strength of global synchronization in scale-free ECoG during the eyes-open, eyes-closed, and sleep states. Specifically, a global correlation map was obtained by calculating the cross-correlation between the BBP fluctuation at every channel and the global BBP fluctuation obtained by averaging across all channels. Note that we referred to the latter as the “global signal” because the global resting-state fMRI activity is defined in the same way. The resulting channel-wise correlation coefficients (converted to the z-score) were further averaged across channels to yield a single measure of the strength of global synchronization in scale-free activity for each session and each arousal state. We then compared the differences in the strength of global synchronization between the sleep and eyes-closed states and between the eyes-closed and eyes-open states. The significance of such differences was assessed by paired t test with p < 0.05.
Cross-correlation between scale-free EEG and global fMRI.
We used concurrently acquired fMRI and EEG data to analyze the coupling between scale-free electrophysiology and global fMRI. From each EEG sensor, the scale-free and the oscillatory spectrograms were separated using the aforementioned IRASA method. Averaging the spectrograms from all sensors yielded the “global” scale-free spectrogram and the oscillatory spectrogram. From these global spectrograms, two power fluctuations were obtained for any specific frequency point: one from the scale-free component and the other from the oscillatory component. Both power fluctuations were resampled to the time points of the fMRI signal. The global fMRI signal was obtained by averaging the signals from all brain voxels. The cross-correlation between the global fMRI and the global power fluctuation was evaluated at every frequency for either the scale-free or oscillatory component of EEG with varying time lags from −30 to 30 s. The resulting correlations were first calculated for each subject and then averaged across subjects.
For the scale-free component of EEG, the frequency-specific cross-correlation functions were further averaged across all frequencies to yield a transfer function that presumably linked the scale-free EEG to the global fMRI. For comparison, another transfer function was also estimated in an attempt to link the oscillatory EEG to the global fMRI. The difference between these two transfer functions was tested for significance for each time point between −30 and 30 s using paired t test across 19 subjects (DOF = 18, p < 0.001 without correction). With this statistical test, we aimed to address whether a frequency-independent coupling between EEG and fMRI is unique for the scale-free component of EEG, but not for its oscillatory component.
In addition, we calculated the cross-correlation between the global fMRI signal and the broadband power fluctuation of the scale-free component at every EEG channel. Such correlations were calculated at 2 time lags (5 s and 12.5 s) by which fMRI was set to be delayed from EEG. The significance of such correlations (converted to the z-score) was evaluated channel by channel using one-sample t test across 19 subjects (DOF = 18, p < 0.005, Bonferroni corrected for the total number of channels). For comparison, we also evaluated the cross-correlation between the global fMRI and the BLP fluctuation.
To further examine whether scale-free EEG contributes to the fMRI signal globally, we evaluated the cross-correlation between the global broadband power fluctuation of scale-free EEG and the fMRI signal of every voxel with two time lags (5 and 12.5 s). Statistical significance of the cross-correlation (converted to the z-score) was tested for every voxel using the one-sample t test (DOF = 18, p < 0.005 without correction). For comparison, cross-correlation was also evaluated between the average alpha-band power fluctuation (obtained by averaging the alpha BLP fluctuations across channels) and the fMRI signal. We also compared the percentage of voxels where the fMRI signals were significantly correlated with the average BBP signal with the average alpha-BLP signal and tested the significance of this difference using nonparametric Wilcoxon signed-rank test with p < 0.01.
Results
We used a recently developed IRASA method (Wen and Liu, 2016) to separate the broadband scale-free and band-limited oscillatory components from the power spectrum of electrophysiological activity observed with macaque ECoG (n = 2), and human MEG (n = 3) or EEG (n = 19) acquired simultaneously during fMRI.
Power-law distribution of arrhythmic electrophysiology
Figure 1A shows the separated scale-free and oscillatory power spectra for an example ECoG signal from a macaque brain during eyes-closed wakefulness. The oscillatory component revealed prominent delta (0–4 Hz) and alpha (8–13 Hz) rhythms. The scale-free component showed a broadband power distribution that could be modeled by a piecewise power-law function. In the log–log scale, this function was shown as two descending lines intersecting at a “knee” frequency of 15 ± 4 Hz across channels, modalities, and arousal states. Therefore, we set out to characterize the broadband scale-free component separately for the LF (1–15 Hz) and HF (15–100 Hz) ranges using BBP, which was the average log power quantified based on the power-law fit of the broadband spectral distribution in the corresponding frequency range (Miller et al., 2009; Wen and Liu, 2016).
Separation of oscillatory and scale-free ECoG. A, Scale-free and oscillatory components of macaque ECoG signals were separated into broadband and narrow-band spectral distributions in the frequency domain. B, From the oscillatory spectrogram, band-limited power fluctuations were extracted from a specific frequency of interest (e.g., alpha). From the scale-free spectrogram, broadband power fluctuations were extracted from either the LF range (1–15 Hz) or the HF range (15–100 Hz).
Slow-power fluctuation of scale-free neural activity
To examine both narrow-band and broadband spectral changes over time, we applied IRASA to the ECoG signal in a 3 s sliding window with a 1 s increment, yielding the separated time–frequency representations (i.e., spectrograms) of the oscillatory and scale-free components (Fig. 1B). For the oscillatory component, the spectral change was frequency specific, giving rise to BLP signals with characteristic carrier frequencies (e.g., alpha). For the scale-free component, the powers at individual frequencies varied in a coherent way that followed the power-law relationship. The fluctuation level of BBP was 11.9 ± 3% (or 12.1 ± 4%) for the LF (or HF) range during sleep, whereas it was 7.94 ± 1.7% (or 9.1 ± 3.6%) during eyes-closed wakefulness. Significant differences (p < 0.01, two-sample F test with Bonferroni correction) between these two states in terms of their LF (or HF) BBP fluctuation levels were observed for 85 ± 4% (or 73 ± 13%) channels. Because both BLP and BBP fluctuated with similar time scales (mostly >10 s) as spontaneous fMRI signals, the scale-free and oscillatory component of electrophysiology were both possible neural correlates to resting-state fMRI.
Scale-free ECoG and MEG are globally correlated
We set out to evaluate the interregional temporal correlation in oscillatory or scale-free spectral fluctuations for various arousal states, including eyes-open and eyes-closed wakefulness and sleep. With 128-channel macaque ECoG recordings from the lateral cortical surface, we first conducted a seed-based correlation analysis similar to the one used widely for resting-state fMRI. For a seed location in the motor cortex, as an example, the alpha-BLP signals were correlated within the local area or network, whereas the LF-BBP signals were globally synchronized (Fig. 2A). For LF BBP, 87 ± 9% of channels were significantly correlated (p < 0.001, Bonferroni corrected) with the selected seed channel, significantly greater than 19 ± 15% for alpha BLP (p = 0.0279, paired Wilcoxon signed-rank test). For any seed channel, ∼80% of all pairs of channels were significantly correlated (p < 0.001, Bonferroni corrected) in BBP during wakefulness and further increased to >90% during sleep, being significantly higher than those for BLP (p < 0.01 during wakefulness, p < 0.05 during sleep, paired Wilcoxon signed-rank test; Fig. 2B). The pairwise correlations in BBP were significantly less inhomogeneous than those in BLP (p < 0.001, paired t test; Fig. 2B). In addition, the correlations in BBP were significantly more global and nonspecific during sleep than during eyes-open or eyes-closed wakefulness (p < 0.05, paired Wilcoxon signed-rank test for percentage and t test for inhomogeneity; Fig. 2B).
Broadband power correlations for scale-free ECoG. A, Seed-based correlation patterns of the LF BBP and the alpha BLP in eyes-closed wakefulness for macaque C. B, Percentage of significantly correlated pairs of channels (left) and the inhomogeneity of cross-correlations between different channels. The asterisk bracket indicates the statistically significant difference between samples (*p < 0.05, **p < 0.01, ***p < 0.001, paired Wilcoxon signed-rank test for percentage and paired t test for inhomogeneity. C, Pattern of correlation in the LF BBP between every channel and the average across all channels. D, Degree of global correlation across three arousal states. Each line represents one session. Relative to sleep, the global correlation in eyes-closed wakefulness was significantly lower (p = 0.014, paired t test). The global correlation in eyes-open wakefulness was further lower (p = 0.018, paired t test) than that in eyes-closed wakefulness.
To further assess the level of global synchrony independently of the seed location, we calculated the cross-correlation in the LF-BBP signal between every ECoG channel and the average across all channels. Such correlations to the global signal were generally high and widely distributed. Although the pattern of global synchrony was preserved across different arousal states (Fig. 2C), the level of global synchrony varied across states, being stronger during sleep than during wakefulness with eyes closed or open (Fig. 2D). Compared with sleep, the global correlation in eyes-closed wakefulness was lower by Δr = −0.137 (p = 0.014, paired t test); relative to eyes-closed wakefulness, the global correlation in eyes-open wakefulness was further lower by Δr = −0.048 (p = 0.018, paired t test). Similar findings were also observed with the HF-BBP. The HF-BBP and LF-BBP were also found to covary at many ECoG channels over nearly the entire lateral surface of the cortex in eyes-open, eyes-closed, and sleep states (Fig. 3A,B), suggesting a common source of modulation affecting both frequency ranges.
Cofluctuations of LF-BBP and HF-BBP. A, BBP fluctuations in the LF and HF ranges were correlated by r = 0.47 for one typical ECoG channel. B, Correlations between LF-BBP and HF-BBP for all 128 sensors were displayed as a map (left). Right, Percentage of channels showing a significant correlation between LF-BBP and HF-BBP fluctuations for various arousal states.
We repeated the above analysis for 274-channel whole-head MEG recordings when human subjects naturally fell asleep from eyes-open and eyes-closed wakefulness. As found with macaque ECoG, the LF and HF BBP were found to covary at 98% of channels in sleep state. The LF (or HF) BBP of scale-free MEG activity not only varied across arousal states, but also fluctuated within each state to a stronger degree during sleep than during wakefulness (Fig. 4A). The fluctuation level of LF-BBP or HF-BBP was significantly higher (p < 0.01, two-sample F test with Bonferroni correction) in the sleep state (9 ± 0.9% or 7 ± 1%) than in the eyes-open state (6.8 ± 0.8% or 4.7 ± 0.7%) or eyes-closed awake state (7 ± 1% or 5 ± 1.3%) for 98% of channels. In addition, the correlation structure and pattern tended to be global for the power fluctuation of the scale-free component but regional for the oscillatory component (Fig. 4B). During sleep, 97 ± 1% (or 98 ± 1.6%) of the total pairs of different MEG channels were significantly correlated in LF-BBP (or HF-BBP), whereas it was 41 ± 9% for BLP. The pairwise correlations in LF-BPP (or HF-BBP) were also more nonspecific, showing significantly lower levels of inhomogeneity than that for BLP (p = 0.0024 for LF-BBP vs BLP or p = 0.004 for HF-BBP vs BLP, paired t test).
Power fluctuations and correlations of oscillatory and scale-free MEG. A, Alpha BLP, and LF and HF BBP fluctuations (right) were extracted from the spectrograms of the oscillatory and scale-free components (left), respectively. B, Patterns of cross-channel correlations in the alpha BLP or the LF BBP. Each row is a seed-based correlation map displayed as a 2D surface topography.
In interpreting the above results, it was possible that the observed correlational differences might reflect the difference in signal-to-noise ratio (SNR) between different signals or states instead of a genuine signal relationship. Although this possibility was difficult to rule out, our observations could not be entirely due to the difference in SNR. We compared the amplitudes of BBP and BLP during every session and every state for both ECoG and MEG. The amplitudes of HF-BBP were significantly lower than those of alpha-BLP (p < 10−4, paired t test), implying that SNR was lower for HF-BBP than for alpha-BLP. Such a difference in SNR was opposite to the observed difference in correlation because HF-BBP was more globally and nonspecifically correlated across channels than was alpha-BLP. In addition, the amplitude of HF-BBP showed no significant difference across awake and sleep states (p = 0.86, paired t test), suggesting that different states had comparable SNRs that could not account for the aforementioned observation that the degree of global correlation was stronger during sleep than wakefulness.
Together, these results from macaque ECoG and human MEG data suggest that the temporal fluctuation of scale-free electrophysiological activity is globally and nonspecifically synchronized to a varying degree that depends on the arousal level or state.
Scale-free EEG is coupled with global resting-state fMRI
Given the above findings, we further hypothesized that the scale-free electrophysiological signal reflected the neural origin of the global fMRI signal. To test this hypothesis, we evaluated the coupling between concurrently acquired fMRI and EEG data in the eyes-closed resting state. As in the ECoG and MEG analyses, we extracted separate spectrograms for the oscillatory and scale-free neuroelectric components from EEG and averaged them across channels (see Fig. 5A for an example). For either component, we calculated the cross-correlation between the power fluctuation at every frequency and the global fMRI signal obtained by averaging the time series of all brain voxels. The resulting correlation coefficients, described as a function of EEG frequency (up to 100 Hz) and time lags (−30 to 30 s) between fMRI and EEG, were low and variable across frequencies for the oscillatory component, but significantly higher and consistent across frequencies for the scale-free component (Fig. 5B). This result indicates the tight coupling between scale-free EEG and global fMRI. This coupling could be further represented by a linear transfer function obtained by averaging the time-shifted cross-correlations across all frequencies. This transfer function implies that an impulse change in scale-free EEG is expected to cause a delayed global fMRI response with a narrow positive peak at ∼5 s and a broad negative peak at ∼12.5 s (Fig. 5C). This transfer function generally agrees with the hemodynamic response function (HRF) used to model the neurovascular coupling and suggests a causal relationship between scale-free EEG and global fMRI. In addition, we also evaluated the cross-correlation between the global fMRI signal and the BBP or BLP signal at every EEG channel. We found that the global fMRI was globally and significantly correlated with EEG BBP, but not with BLP, when fMRI was delayed from EEG by 5 s or 12.5 s (p < 0.005, Bonferroni correction for the number of channels; Fig. 6A,B). Note that the transfer function that linked scale-free EEG to global fMRI, although qualitatively similar to the HRF, had a deeper undershoot peaked at 12.5 s (Fig. 5) and the integral of the transfer function was slightly negative. The implication of this observation is unclear and the underlying mechanism is entirely mysterious.
Direct coupling between scale-free EEG and global fMRI. A, Oscillatory and scale-free spectrograms observed with EEG. B, Cross-correlation between the global fMRI and the spectral power fluctuation at every frequency with varying time delay for the oscillatory (top) and scale-free components (bottom). C, Cross-correlations were averaged across frequencies to yield a linear transfer function to relate oscillatory (blue) or scale-free (red) EEG fluctuation to global fMRI activity. The asterisk bracket indicates the statistically significant difference between the two correlation functions within the specified ranges of time delay (*p < 0.001, paired t test).
Correlation between fMRI signals and the power fluctuations of scale-free and oscillatory EEG. A, B, Cross-correlations between LF-BBP, HF-BBP, alpha-BLP at every channel and the global fMRI signal as fMRI delayed from EEG by 5 s (A) or 12.5 s (B). Channels with significant or nonsignificant correlation are marked as dark versus gray dots. C, D, Cross-correlations between the global LF-BBP, HF-BBP, alpha-BLP, and the fMRI signal at every voxel as fMRI delayed from EEG by 5 s (C) or 12.5 s (D). Colored bars indicate the correlation coefficient after the Fisher's transformation.
In addition, we also evaluated the cross-correlation between the fMRI signal at every voxel and the global BBP fluctuation averaged across all EEG channels. We found that the LF and HF BBP fluctuations were significantly and positively (or negatively) correlated with the fMRI signals when fMRI was delayed from EEG by 5 s (or 12.5 s) for 69% (or 64%) and 65% (or 70%) of the cortical gray matter, respectively (Fig. 6C,D). In contrast, the correlation pattern between alpha-BLP and fMRI at a 5 s delay was confined to specific brain regions, being positive in the thalamus but negative in the visual cortex. Significant fMRI-BLP correlations were found at ∼14% of the cortical voxels, which was significantly lower than those from BBP (p = 0.0034 for LF-BBP vs BLP, or p = 0.0016 for HF-BBP vs BLP, paired Wilcoxon signed-rank test; Fig. 6C). The above findings support the notion that scale-free EEG is the neuroelectric correlate of global fMRI fluctuation in the resting state.
Discussion
We investigated the neural origins and correlates of resting-state fMRI. Our results demonstrate the following: (1) scale-free and oscillatory components could be separated from neurophysiological signals; (2) broadband power fluctuations of scale-free ECoG and MEG were globally correlated in various arousal states; (3) broadband fluctuation of scale-free EEG was directly coupled with the concurrent global fMRI signal; and (4) these spatiotemporal properties of scale-free neural activity in relation to resting-state fMRI were notably different from those of band-limited neural oscillations, of which the power fluctuations were correlated within specific regions or networks. Therefore, scale-free and oscillatory neural processes both contribute to spontaneous fMRI signals, but account selectively for global and modular spatial patterns of functional connectivity, respectively.
Multitude neural correlates to resting-state fMRI
Resting-state fMRI has a complex and elusive basis related to a multitude of neuronal and non-neuronal processes (Leopold and Maier, 2012). The fMRI signal has been found to correlate to the power fluctuations of various frequency components of electrophysiology both during tasks (Mukamel et al., 2005; Niessing et al., 2005; Goense and Logothetis, 2008; Scheeringa et al., 2011) and at rest (Mantini et al., 2007; Shmuel and Leopold, 2008; de Munck et al., 2009; Schölvinck et al., 2010; Liu et al., 2014; Hipp and Siegel, 2015). Therefore, fMRI is shaped collectively by all frequency components instead of indicating a single frequency component.
Different frequency components may not always arise from distinct neural processes. In fact, neural activity exhibits a variable mixture of arrhythmic and rhythmic dynamics (He et al., 2010; Buzsáki et al., 2012). These two types of dynamics differ in their spectral profiles: the arrhythmic activity is broadband and scale free, whereas the rhythmic activity is narrow-band and scale specific. Without dissociating them, one has to limit themselves to time periods or frequency ranges with either rhythmic or arrhythmic activity, but not both, to study their individual dynamics while being less concerned about their mutual interference.
In contrast, we advocate a two-step strategy. First, neural signals are separated into a broadband scale-free component and a set of narrow-band oscillatory components. The IRASA method used here (Wen and Liu, 2016) and other methods (Yamamoto and Hughson, 1991; He et al., 2010; Miller et al., 2014) provide the tools with which to dissociate oscillations at distinct frequencies from frequency-independent broadband activity and to isolate arrhythmic activity from coexisting oscillations. Second, we can characterize the separated scale-free or oscillatory components and then compare them with the fMRI signal to uncover their common features and direct relationships. However, IRASA is unable to separate scale-free and oscillatory components in the time domain or to uncover their phase–phase or phase–amplitude relationship.
Scale-free neural activity underlies global resting-state fMRI activity
Scale-free neurophysiological signals can be well described by a power-law function of frequency (Bédard et al., 2006; Manning et al., 2009; He et al., 2010; Miller et al., 2014). Individual frequencies are mutually dependent as constrained by the power-law distribution. When the power-law function varies over time, the powers at distinct frequencies cofluctuate accordingly. Therefore, scale-free neural activity is uniquely suited to explain previous findings that correlations between fMRI signal and LFP power were consistent across frequencies within a broad range (Goense and Logothesis, 2008; Schölvinck et al., 2010). For its broadband nature, it is also more appropriate to characterize scale-free activity with its broadband power than the (likely collinear) band-limited powers at individual frequencies.
We showed here that scale-free neural signals exhibited spontaneous power fluctuations that were globally synchronized and directly coupled with global resting-state fMRI activity. These findings extended the previous finding that widespread resting-state fMRI activity was correlated with the LFP power fluctuation at a single cortical site (Schölvinck et al., 2010). First, we further specified the neural basis of global fMRI activity to scale-free neural activity. The broadband distribution of scale-free activity could well explain the consistent correlations between the global fMRI signal and the band-limited powers of LFP across frequencies in the gamma band (>40 Hz) in the absence of oscillations (Schölvinck et al., 2010). Such broadband correlations could be generalized to the entire spectral range (Fig. 5) after using IRASA to exclude oscillations. Second, distributed neural recordings with ECoG and MEG offered additional evidence that scale-free activity was indeed globally correlated (Figs. 2, 4). It further confirms its role as the neural basis of global resting-state fMRI activity. Last, neural activity underlying global fMRI should be widespread in itself and thus be observable across large spatial scales. We showed that the neural basis of global resting-state fMRI activity also manifested itself as synchronized population activities measured with ECoG, MEG, and EEG for both macaques and humans. It further calls for serious caution in removing the global signal in resting-state fMRI analysis (Murphy et al., 2009; Saad et al., 2012).
State dependence of widespread scale-free activity
The fluctuation and correlation in the broadband power of scale-free ECoG was more intensive and global during sleep than wakefulness (Fig. 2), suggesting the dependence of scale-free activity on arousal states. Spontaneous fMRI signals were similarly different in sleep versus wakeful states (Fukunaga et al., 2006; Horovitz et al., 2008). There is thus likely a vigilance-related mechanism that supports the widespread and nonspecific fluctuation in both scale-free electrophysiology and fMRI.
Indeed, the reticular activation system (RAS) governs the awareness, arousal, and vigilance level and innervates the entire cortex through diffuse neuromodulation pathways (Magoun, 1952). These pathways comprise a structural network to ascend time-variant neurochemical (e.g., cholinergic, Mesulam, 1995), common inputs to drive the global cortical fluctuation while regulating cortical state (Harris and Thiele, 2011) and sleep–wake transition (Saper et al., 2001). Stimulating RAS disrupts cortical oscillations and induces arrhythmic activity, indicating the role of RAS in modulating scale-free electrophysiology (Moruzzi and Magoun, 1949; Munk et al., 1996). Other recent studies also emphasize the vigilance-related physiological processes or events as the source of global activity fluctuations (Wong et al., 2013; Liu et al., 2015). The specific relationship among RAS, scale-free activity, and vigilance remains speculative, and needs to be elucidated in future studies.
Neural and hemodynamic activity across multiple scales
LFP, ECoG, MEG, and EEG result from the sum of population activity within increasingly larger spatial scales. Their precise relations to underlying neuronal activity are complex (Buzsáki et al., 2012). fMRI measures hemodynamic activity in a spatial scale closer to that of LFP but considerably different from those of ECoG/EEG/MEG. The temporal relationship between neural and hemodynamic signals across different spatial scales is not straightforward to interpret.
Spiking activity is not directly observable with far-field recordings (e.g., ECoG, MEG, and EEG), whereas synaptic activity is observable with LFP and far-field ECoG and MEG/EEG if they are aligned spatially (Buzsáki et al., 2012). Importantly, spiking activity has been shown to correlate with the high-gamma or broadband power of synaptic activity observed with LFP (Ray and Maunsell, 2011), ECoG (Manning et al., 2009; Miller et al., 2009, 2014) and EEG (Whittingstall and Logothesis, 2009). It is plausible that presynaptic spiking activity arrives in synchrony to a population of spatially aligned neurons such that the vector sum of postsynaptic currents is strong enough to produce macroscopic signals. When neuronal spikes are temporally irregular and arrhythmic, population-averaged presynaptic spiking activity and postsynaptic current activity should both be (time) scale free and exhibit a power-law spectral distribution. Therefore, the broadband power of population synaptic activity is expected to correlate to spiking activity averaged in the period when the power spectrum is evaluated. Because the fMRI signal is coupled to population synaptic activity (Logothetis et al., 2001), this speculative scenario further implies that the fMRI signal is correlated with both broadband synaptic activity and arrhythmic spiking activity at the level of cortical populations.
Neural oscillations versus resting-state fMRI
Oscillations have been common targets for neurophysiological studies (Buzsáki and Draguhn, 2004; Buzsáki, 2006). Oscillations may contribute to network fMRI activity in an informative manner (Hipp et al., 2012; Siegel et al., 2012), but be obscured or confounded by broadband neural activity that also correlates to fMRI. This confound is of particular concern for the low-frequency range (<30 Hz), in which oscillations are most notable (Manning et al., 2009; Miller et al., 2014). This may explain previous contradicting findings that changes in the delta, alpha, or beta power could be correlated with fMRI both negatively and positively across various conditions and regions (Mukamel et al., 2005; Goense and Logothetis, 2008; de Munck et al., 2009).
After removing the scale-free component from neural signals, the remaining oscillatory component was truly band limited and timescale specific. Unlike scale-free activity, the power fluctuations in (alpha) oscillatory activity were correlated within specific regions or networks and the correlations between the power of (alpha) oscillatory activity and fMRI signals was confined to specific thalamic and cortical sites. These results led us to deduce that scale-free and oscillatory neural activities, while both contributing to hemodynamic fluctuations, account for distinct spatial patterns of correlations in resting-state fMRI signals. Specifically, globally and modularly orchestrated RSNs report scale-free and oscillatory electrophysiological processes and interactions, respectively.
Footnotes
This work was supported in part by the National Institutes of Health (Grant R01MH104402) and the Purdue Research Foundation. We thank Dr. Naotaka Fujii for publicizing the macaque ECoG data and Drs. Jeff Duyn, Masaki Fukunaga, and Jacco de Zwart for assistance in collecting MEG and fMRI-EEG data.
- Correspondence should be addressed to Zhongming Liu, PhD, Assistant Professor of Biomedical Engineering, Assistant Professor of Electrical and Computer Engineering, College of Engineering, Purdue University, 206 S. Martin Jischke Dr., West Lafayette, IN 47907. zmliu{at}purdue.edu