The temporal structures and functional significance of scale-free brain activity

Power Spectra of Electrical Field Potentials of the Brain

We collected spontaneous electrocorticography (ECoG, i.e., invasive EEG) data from five patients (Patients #1 - #5, Table S1) undergoing surgical treatment for drug-resistant epilepsy, in both quiet wakefulness and slow-wave sleep (SWS, sleep stages 3/4). The length of data collected in each arousal state ranged from 12 to 83 min. In Patient #4, 30 min of rapid-eye-movement (REM) sleep data were also collected. Electrode locations documented by plain X-ray pictures are shown in Fig. 1F. For the details of clinical and data collection information see Table S1.

The power spectra of ECoG signals, plotted in log-log coordinates (Fig. 1 A-E, top graphs), roughly followed a straight line, with local peaks corresponding to well-known brain oscillations rising above this line. These oscillations included the slow oscillation (~0.8 Hz) and sleep spindles (~12 Hz) during SWS, the theta (~7-8 Hz), alpha (~10 Hz) and beta (~20 Hz) oscillations in the awake state. Multiple brain oscillations notwithstanding, the vast majority of power in the ECoG signals is accounted for by the power-law distribution representing arrhythmic brain activity.

To estimate the power-law exponent, we first adopted the procedure developed by Yamamoto and Hughson (1991, 1993), coarse-graining spectral analysis (CGSA), to separate the harmonic/oscillatory from the scale-free/fractal components of the ECoG power spectrum. This method takes advantage of the self-affinity property of a scale-free time series, i.e., the statistical distribution of the data remains the same when sampled at different scales (Mandelbrot and Van Ness, 1968), and the fact that this is not true for harmonic time series (for details see Supplementary Experimental Procedures). The fractal and harmonic components of the ECoG power spectrum obtained by the CGSA method for each patient, each arousal state are presented in Fig. 1 A-E (middle graphs). In SWS, harmonic components included peaks at 0.5-1 Hz and 12-14 Hz, representing the slow oscillation (also called “up-and-down states”, UDS) and sleep spindles respectively. In the waking state, the harmonic components included various peaks in the delta, theta, alpha, beta and gamma frequency ranges. In the REM sleep data from Patient #4, the harmonic component included peaks in the theta frequency range and around 30 Hz.

Interestingly, after reducing oscillatory activity, the remaining fractal components of the power spectrum still had a “shoulder” in the frequency range of 0.1-1 Hz. Whether this reflected a true property of arrhythmic brain activity, or alternatively, the remaining presence of the UDS that was unsuccessfully removed by the CGSA method awaits future investigation². Accordingly, we fit the low-frequency end (<0.1 Hz) and the higher-frequency range (1-100 Hz) respectively with a power-law function: P ∝ f^-β, avoiding the “shoulder”. The values of the resultant β exponents are listed in Fig. 1A-E. Averaged across all patients, the low (<0.1 Hz) and high (1-100 Hz) frequency ranges had exponent of 2.2 and 2.44 respectively in the awake state, and 1.6 and 2.87 respectively in SWS. The effect of arousal state on power-law exponent was non-significant (P > 0.5). However, the interaction of arousal state × frequency range was highly significant (P < 0.004); whether this reflected unreduced UDS in the data or a true change in arrhythmic brain activity across waking and SWS requires future investigation. Notably, the exponents we obtained correspond very well with previous reports on the power-law exponent of LFP or ECoG activity (Freeman and Zhai, 2009; Miller et al., 2009a; Milstein et al., 2009).

The Widespread Presence of Nested Frequencies

To investigate temporal structures within spontaneous ECoG signals, we examined nested-frequency patterns in each electrode. For a pair of frequencies, instantaneous phase and amplitude were extracted for the lower and higher frequency respectively. The lower-frequency phase at all samples was sorted into 20 bins, and the concurrent higher-frequency amplitude was averaged within each bin. The resultant curve shows the dependency of the higher-frequency amplitude on lower-frequency phase, and its deviation from a uniform distribution was evaluated using an inverted entropy measure to yield a modulation index (MI) (Tort et al., 2008). The MI was compared with shuffled data to obtain a MI Z-score indexing the strength of cross-frequency phase-amplitude coupling. For details see Supplementary Experimental Procedures.

Cross-frequency phase-amplitude coupling was investigated in a 2-D frequency space: phase was extracted from each 1-Hz-width bin centered at 1, 2 … 20 Hz, and amplitude was extracted from each 5-Hz-width bin centered at 5, 10 … 200 Hz. For each frequency pair, the percentage of electrodes with a significant MI Z-score (P < 0.05 after Bonferroni correction) is plotted as color in Fig. 1 A-E (bottom graphs). During wakefulness, extensive nested frequencies are present in every patient in the majority of electrodes; the extent of nested frequencies decreases during SWS but remains highly robust across wide frequency ranges in Patients #3, #4 & #5. REM sleep data from Patient #4 also shows extensive nested frequencies across the entire frequency spectrum. Such widespread presence of nested frequencies extends beyond the frequency boundaries of rhythmic brain oscillations (see insets in the middle graphs of Fig. 1A-E) and is difficult to reconcile with previous notions conceiving nested frequencies only within the framework of oscillatory brain activities. Instead, it appears that to explain its presence, arrhythmic brain activity must be invoked.

Power-Law Distribution Is Not an Artifact of Averaging

Brain oscillations and “1/f noise” have so far remained two largely separate fields. An important, unresolved question concerns whether they index distinct physiological processes with different underlying mechanisms (e.g., Bullock et al., 2003; Freeman and Zhai, 2009), or alternatively, if “1/f noise” is the result of averaging across time of many different, transient and recurrent oscillations at different frequencies and with different amplitudes (e.g., Buzsaki, 2006).

To investigate this issue, we randomly picked three electrodes (from Patient #3), two with and one without rhythmic oscillations. As shown by the power spectra averaged across the entire awake record (Fig. 2A left), electrode #33 contains oscillations at ~1.5 Hz and ~20 Hz, electrode #43 contains oscillations at 7-8 Hz and ~20 Hz, and electrode #64 contains no periodic oscillations but only arrhythmic, scale-free activity. The above impression from the power spectra is corroborated by visual inspection of the raw ECoG records. Two randomly selected 20-second segments of raw data, which were separated in time by >1 hour are shown in Fig. 2C. The power spectra of these two short data segments (Fig. 2A middle and right) recapitulate the power spectra averaged over the entire 83-min record. The observation that the overall shape of the power spectrum – its power-law distribution and the location of oscillatory peaks – appears to be stable over time was obtained in many electrodes.

These results suggest that the power-law distribution of the ECoG power spectrum does not seem to be an artifact of averaging over time of many independent, periodic oscillations, but rather appears to be the direct result of the presence of arrhythmic, scale-free brain activity. One remaining possibility is that there are hidden rhythmic activities at finer spatial scales than the resolution of ECoG electrodes. However, arrhythmic, scale-free activity is present at every spatial scale in the brain. The power spectrum of spontaneous LFP activity also follows a power-law distribution (Leopold et al., 2003; Milstein et al., 2009) with a power-law exponent close to that of ECoG activity. Further, cortical pyramidal neurons’ spikes are usually nonrhythmic (Ermentrout et al., 2008; Faisal et al., 2008; Freeman and Zhai, 2009; Koch, 1997; Thivierge and Cisek, 2008; but see Maimon and Assad, 2009) and “1/f noise” has been observed in neuronal spike trains (Gisiger, 2001; Gruneis et al., 1989; Takahashi et al., 2004; Yamamoto, 1991) and neurotransmitter release (Lowen et al., 1997).

Stability of Nested Frequencies within Arrhythmic Brain Activity

Next we investigated whether electrodes with only scale-free brain activity but no rhythmic oscillations also contained nested frequencies. Nested frequencies in the above three electrodes were computed in the same 2-D frequency space as shown in Fig. 1A-E (bottom graphs). Significant cross-frequency coupling across wide frequency ranges were found in every electrode (Fig. 2B), including electrode #64, which contained no discernible periodic oscillations.

To examine the stability of nested-frequency patterns within arrhythmic brain activity, we plotted the higher-frequency amplitude averaged at different phases of the lower frequency for selected frequency pairs in electrode #64. Amplitude of 5-Hz-width bands centered at 25, 50, 100, 150, 200 Hz was each averaged at different phases of lower-frequency bands (1-Hz-width) centered at 1, 6, 11, 16 Hz and plotted in Fig. 2D. This analysis was repeated for the entire awake record (Fig. 2D, left) and the two 20-sec segments shown in Fig. 2C (Fig. 2D, middle and right). The nested-frequency patterns from the two short data segments are each similar to the pattern from the entire record, suggesting that the nested frequencies present within arrhythmic brain activity appear to be stable over time. Similar results from a few other representative electrodes in each patient are shown in Fig. S1. Importantly, the pattern of nested frequencies varied across electrodes, even though all of them contained mainly arrhythmic brain activity (Fig. S1).

Nested-Frequency Patterns across Brain Surface

To characterize nested-frequency patterns across electrodes, for the same frequency pairs as in Fig. 2D, we obtained the MI Z-score, which indexes the strength of cross-frequency phase-amplitude coupling, for each of the 259 electrodes from all five patients. The percentage of electrodes with significant cross-frequency phase-amplitude coupling for each frequency pair is listed in Table 1 (the bottom values in each cell). Averaged across the 20 frequency pairs, 85% and 72% electrodes had a significant MI Z-score in wakefulness and SWS respectively.

For each electrode and frequency pair we also determined the preferred phase of the lower frequency, which denotes the phase of the lower frequency fluctuation at which the higher frequency has the largest amplitude. Plotting the MI Z-score against preferred phase on an electrode-by-electrode basis for each patient and frequency pair, we found that the preferred phase tends to cluster around 0 and ±π, that is, the peak and trough of the lower-frequency fluctuation respectively (Figs. 3A & S2). For each frequency pair, the percentage of electrodes (across all five patients) with preferred lower-frequency phase around 0 or ±π, and significant cross-frequency coupling, are listed in Table 1.

The above bimodal distribution of preferred phase raises the question of whether the electrodes with preferred phase around 0 and those with preferred phase around ±π situated in different cortical regions. Thus, for each frequency pair, we plotted the preferred phase of the lower frequency as color on a 2-D representation of the electrode grid (Patient #3 in Fig. 3B, Patients #1, 4&5 in Fig. S3). Indeed, electrodes preferring phase 0 and those preferring phase ±π formed largely separate clusters, which to a rough degree followed cortical anatomy such as the location of Sylvian fissure and the central sulcus (marked in Fig. 3B). Moreover, these spatial patterns were relatively stable across arousal states (significant spatial correlations across wakefulness and SWS are marked by squares in Figs. 3B and S3).

It should be noted that the results of this analysis were contributed by both arrhythmic brain activity and periodic brain oscillations, because, unlike in Fig. 2, spontaneous ECoG signals from all electrodes were used. When periodic brain oscillations and arrhythmic brain activity are already mixed in collected data, as is the case for many electrodes in our recordings, there is no a posteriori mathematical means to separate them while preserving their respective biologically-meaningful phase. The influence of brain oscillations on the results in Figs. 3 & S2 can be most clearly seen in 1-Hz-phase modulation of higher frequency amplitudes (the leftmost column in each sub-panel, also see Table 1). During wakefulness, the preferred 1-Hz phase is a bimodal distribution, which clustered around 0 and ±π; however, during SWS it becomes a unimodal distribution with most electrodes having their preferred phase around phase 0. This is in full accordance with known neurophysiology: During SWS, the slow oscillation (i.e., UDS) at around 1 Hz becomes very prominent, which modulates higher frequencies strongly. The surface positive (phase ~0) ECoG activity at ~1 Hz corresponds to the intracellular “up-state”, during which higher-frequency activities are dramatically increased (Vyazovskiy et al., 2009).

Task Modulation of Scale-Free Brain Activity

We tested whether scale-free brain activity has any functional significance by recording ECoG signals during both quiet wakefulness and task performance in Patients #3 & #5 and three additional subjects (#6 - #8, see Table S1). The task included a cued button press condition and a self-paced button press condition. In the cued condition, a visual cue prompted the subject to press a button as soon as they detected the cue, the inter-trial interval varied randomly between 2 and 20 seconds (mean 6.1 sec). In the self-paced condition, the subject was instructed to press the button at their own pace, separating adjacent button presses by a few seconds but avoiding regular rhythms (which could entrain arrhythmic brain activity into oscillations, see, e.g., Elbert et al., 1991; Lakatos et al., 2008). The button press was performed by either the left or right index finger, each corresponding to either the contralateral or ipsilateral hand, depending on which brain hemisphere the electrodes covered.

We searched blindly for electrodes whose power spectrum was significantly altered by task performance. Significant deviation of power-law exponent during task from that during rest was found in 11 electrodes, including four over the hand representation of primary motor cortex, one over lower motor cortex (tongue/face representation), three over premotor area, one over Broca’s area, one over lateral occipital cortex, and one over lateral temporal lobe (Fig. 4 A-C & Fig. S4). All five electrodes over primary motor cortex were confirmed by cortical stimulation as part of patients’ routine clinical care. Compared to rest, the power-law exponent β of all these electrodes decreased during all four trial types (i.e., the slope became “shallower”), suggesting that scale-free brain activity is modulated by task performance. Interestingly, there was no systematic difference between the four trial types. Consistent with a previous report (Miller et al., 2009a), we did not find any significant change of power-law exponent in the higherfrequency range (10-100 Hz) during task performance.

In electrode #64 from Patient #3 that contained only scale-free brain activity during the spontaneous waking state (see Fig. 2), we found the emergence of an oscillation at ~8 Hz during task performance in all four trial types (Fig. 4D). This oscillatory activity can be seen in both the power spectrum and raw data records. Nonetheless, the nested-frequency patterns did not change (compare Fig. 4D bottom right graph to Fig. 2D, and see Fig. S5). Since our task contained no rhythmic structure, the emergence of this 8-Hz oscillation was not a result of entrainment by external stimuli, but reflected intrinsic properties of the underlying neuronal network.

Variation of Power-Law Exponent across Brain Regions Assessed by fMRI

Since previous studies have demonstrated that the power spectrum of spontaneous fMRI signals also follows a power-law distribution (Bullmore et al., 2001; Fox et al., 2007), we investigated the scale-free dynamics in the fMRI signal and its power-law exponent across brain regions, utilizing the whole-brain coverage provided by fMRI. Resting-state fMRI data was acquired in 17 consenting healthy volunteers in a 3T scanner. A set of 31 regions of interest (ROIs) were obtained from previous studies by the authors or published articles (see Table S2), among which there were 10 pairs of homologous brain regions. Other than three ROIs located outside the neocortex – hippocampus, cerebellum and thalamus – the remaining regions resided in five known brain networks (attention, default, motor, saliency, and visual) (Fox and Raichle, 2007).

The normalized power spectrum (i.e., power spectral density, with total power/variance = 1) was computed for the fMRI signal time course from each ROI for each fMRI run, then averaged across runs and across subjects. Power spectra from homologous brain regions were averaged together, yielding one average power spectrum for each of 21 brain regions (Fig. 5A). The power spectra of spontaneous fMRI signals roughly followed a straight line when plotted in loglog coordinates (Fig. 5A), suggesting that spontaneous fMRI signals are also scale-free. The power-law exponents obtained by a least squares fit are listed in Table S2, with a mean of 0.69 across brain regions, which differed significantly from the power-law exponent of ECoG signals (mean 2.2 for <0.1 Hz range in wakefulness) (P < 10^-6). The origin of this difference between ECoG and fMRI signal power-law exponents warrants future investigation.

The power-law exponent of the fMRI signal was found to vary across brain networks. The visual regions, default network and dorsal anterior cingulate cortex (dACC) had the steepest power spectra, characterized by the largest power-law exponent β. On the other end, the cerebellum and hippocampus had the shallowest power spectra, characterized by the smallest power-law exponent. The motor and attention network regions were intermediate. An ANOVA suggested that the effect of brain network on power-law exponents was highly significant (F_5,15 = 5.05, P = 0.006). By contrast, when the total variance of the fMRI signal (listed in Table S2 for each brain region) instead of the power-law exponent was entered into the ANOVA, no significant network effect was found (P = 0.3). Nonetheless, there was a moderate but significant correlation between fMRI signal variance and power-law exponent across brain regions (Fig. 5C). Of note, the small power-law exponents of the cerebellum and hippocampus reflected a proportionally smaller amount of low-frequency activity as compared to the neocortex, consistent with previous neurophysiological observations and considerations of their respective anatomical structures (Bullock and Basar, 1988; He and Raichle, 2009).

Control Recordings and Analyses

Because fluctuations of resistivity in electronic conducting materials can also exhibit “1/fnoise” (Weissman, 1988), it is important to demonstrate that our data were not contributed by instrument noise. To address this, we conducted a dummy ECoG recording in a standard epilepsy in-patient monitoring room without the presence of a patient. The power spectrum of the dummy recording (duration 46 min) is plotted in log-log coordinates in Fig. S6, together with the average ECoG power spectrum from a patient. The magnitude of the dummy recording is far below that of ECoG signals recorded from the patient. Further, the power spectrum of the dummy recording does not follow a power-law distribution. Fitting a power-law function P ∝ 1/f^β to the low (< 0.1 Hz) and high (1-100 Hz) frequency range nonetheless yielded exponent of 0.6 and 1 respectively, which are significantly different from power-law exponents of brain ECoG signals (P < 0.0001 for both frequency ranges). Moreover, no significant nested frequencies were observed in the dummy recording time series.

Our fMRI data were previously published in Fox et al. (2007). In this study the authors acquired fMRI signals from a water phantom in the same scanner, the power spectrum of which was close to white noise (see Fig. S4 in Fox et al., 2007), which differs significantly from the power spectrum of fMRI signals recorded from the brain (see Fig. 5A).

A recent study showed that microsaccades could produce artifactual broadband activity in the EEG signal (Yuval-Greenberg et al., 2008). Although it was suggested that these concerns do not apply to ECoG recordings (Yuval-Greenberg et al., 2008), we confirmed that the current results were not contributed by muscle artifacts related to eye movements by plotting the strength of cross-frequency coupling across the cortical surface (Fig. S7). All ECoG results presented in this paper were obtained under the common-mean-reference montage, therefore artifacts related to eye movements should be most pronounced in electrodes closest to the eyes (Melloni et al., 2009; Yuval-Greenberg et al., 2008). The results in Fig. S7 confirmed that significant nested frequencies were present across wide cortical regions and did not exhibit an anterior-posterior gradient.

Comparison with Other Scale-Free Dynamics in Nature

Thus far we have shown that arrhythmic, scale-free brain activity contains rich temporal structures and is behaviorally relevant. However, since scale-free dynamics are ubiquitous in nature, one potentially revealing question is whether other scale-free dynamics also contain temporal structures and functional meanings and whether their structures are similar to or different from those observed in brain activity. While it is up to economists and physicists to decide whether scale-free dynamics in market fluctuations contributed to the recent economic crash (Cho, 2009) or if similar dynamics can initiate earthquakes, we merely investigated the potential temporal structures present in these scale-free dynamics and compared them with those in scale-free brain activity.

We conducted nested-frequency analysis on spontaneous earth seismic waves (collected over 4 months) and fluctuations of the Dow-Jones Industrial Average Index (collected over 80 years). Both earth seismic waves and stock market fluctuations followed a power-law distribution in their temporal power spectrum (Fig. 6A). Moreover, the power-law exponent of both data types (1.99 for seismic waves and 1.95 for stock market) was close to that of the brain ECoG activity (mean 2.2 during wakefulness for <0.1 Hz range).

Like arrhythmic brain activity, both earth seismic waves and stock market fluctuations contained extensive nested frequencies (Fig. 6B). However, the exact patterns of nested frequencies in these signals differed from those in brain activity. For example, a bimodal nested-frequency distribution with two preferred phases is seldom seen in the ECoG signals (<10% of all electrodes), but was prominent in earth seismic waves. This bimodal distribution in earth seismic waves might be due to different generators in the earth propagating to the same surface location; by contrast, the rarity of a bimodal nested-frequency pattern in the ECoG signal might imply a more homogeneous source of generators. The observation that nested frequencies are prominent in these natural and economic scale-free dynamics reinforces the conclusion that the extensive nested frequencies observed in ECoG signals were contributed primarily by scale-free brain activity and not by periodic brain oscillations.

Simulations of Scale-Free Dynamics

Finally, to better understand the relationship between nested frequencies and scale-free dynamics, we constructed several simulated scale-free time series and examined whether they contained nested frequencies. All simulated time series were analyzed in a similar manner as the ECoG data for temporal power spectrum and nested frequencies.

First, as a control analysis, we performed nested-frequency analysis on a Gaussian-distributed white-noise time series (Fig. 7A, left). This time series does not contain any nested frequencies (Fig. 7A, right). Second, we constructed a scale-free time series by filtering this white-noise time series in the frequency domain to yield a power-law distribution of P ∝ 1/f^β (β = 1.8) while maintaining the random phase of the white-noise time series. This spectrally-generated scale-free time series with a power-law exponent close to that of brain ECoG activity also did not contain any nested frequencies (Fig. 7B).

Next, because the power-law exponent of the ECoG signals was close to 2, which is typical of a random-walk process (i.e., Brownian motion), we constructed a random-walk model by summing over a white noise time series: x(t) = x(t -1) + ε(t), where ε(t) is the same white noise time series as in Fig. 7A. To our surprise, this random-walk process contained statistically significant nested frequencies (Fig. 7E).

We then added a single parameter to the above random-walk model to change it into a first-order autoregressive (AR-1) process: x(t) = 0.9 x(t-1) + ε(t). The addition of this single parameter eliminated nested-frequency patterns (Fig. 7C). Despite the seemingly small change in the model, an important difference between this AR-1 process and the random-walk process in Fig. 7E is that the random-walk process has long-range temporal dependence (also called “long memory”) whereas the AR-1 process does not.

Next, we constructed a model of an aggregate of AR-1 processes (Fig. 7D). The approach of aggregating a few AR-1 processes to model scale-free dynamics (in a restricted frequency range) has been used extensively in economics (Erland and Greenwood, 2007; Granger, 1980), and has been applied with considerable success to modeling the power-law distribution of psychological data, such as the fluctuations of reaction times in a decision making task (Wagenmakers et al., 2004; Ward, 2002). Contrary to our initial conjecture, such a model did not contain any significant nested frequencies (Fig. 7D). Importantly, although a mixture of AR-1 processes could produce the appearance of “1/f noise”, such dynamics do not have genuine long-range temporal dependence that is present in a random-walk process (Wagenmakers et al., 2004).

It might be counterintuitive that a random-walk process contained structures. One possibility is that this is because of the higher-order statistical regularities present in the arithmetic pseudorandom number generator used for generating the white noise input, which are invisible in most applications but produce interesting patterns under the sensitive nested-frequency analysis and when the modeled process contains long memory. To verify this conjecture, we obtained true random numbers³ from physical source (atmospheric noise, www.random.org). These uniformly-distributed random numbers were transformed into a Gaussian-distributed random number series⁴, and used as the input ε(t) to the random-walk model x(t) = x(t-1) + ε(t). The resultant random-walk time series, with a 1/f² power spectrum, did not have any nested frequencies (Fig. 7F), confirming the foregoing conjecture that the nested frequencies in Fig. 7E were a manifestation of the higher-order statistical regularities present in the pseudorandom number generator.

The above simulations suggest that the combination of even very weak structures in the underlying generative process and long-range temporal dependence are sufficient to produce nested frequencies in scale-free dynamics. However, we caution that the class of linear models constructed here constitutes only a subset of all models that can give rise to scale-free dynamics (for reviews see Bak, 1996; Frank, 2009; Mandelbrot, 1999; Ward and Greenwood, 2007; Weissman, 1988). It would be beyond the scope of the present paper to exhaustively investigate these models, but the current results should encourage future theoretical and modeling work to further explore the temporal structures of scale-free dynamics generated by different mechanisms.

The Interpretation of Nested Frequencies in Scale-Free Dynamics

Whereas the interpretation of nested brain oscillations is relatively straightforward (Jensen and Colgin, 2007; Schroeder and Lakatos, 2009), that of nested frequencies in arrhythmic brain activity is less intuitive. In general, nested frequencies in scale-free dynamics also means that at a particular phase of the lower-frequency fluctuation, such as its trough or peak, the amount of higher-frequency activity increases. For example, the trough of surface-recorded slow cortical potentials (SCPs), which constitute the low-frequency end of the 1/f^β power spectrum (He and Raichle, 2009; Monto et al., 2008), is associated with increased cortical excitability and increased amount of higher-frequency activities (Birbaumer et al., 1990; He and Raichle, 2009; Rockstroh et al., 1989; Vanhatalo et al., 2004). This is in line with comments on “1/f noise” in the physics field suggesting that perturbations at long wavelength cause energy dissipations at all length scales (Bak et al., 1987). Furthermore, the non-random distribution of preferred lowerfrequency phase in ECoG activity, which clustered around phase 0 and ±π, is indicative of specific underlying neurobiological mechanisms.

Nonetheless, as pointed out by Kramer et al. (2008), complex waveforms such as edges or triangular waves could also produce significant results in nested-frequency analysis. While artifacts in EEG/ECoG recordings can certainly create sharp edges, we emphasize that an artifact-free ECoG record (such as electrode #64 in Fig. 2C) is also full of complex waveforms that are non-sinusoidal, indeed waveforms that make up scale-free brain activity. These complex waveforms produced by the brain should not be considered artifacts or noise simply because they are non-oscillatory; their functional meanings deserve and demand future investigation. Our view is similar to that of Bullock et al. (2003), who suggested that “the wide-band, apparently stochastic background activity in the EEG may in part consist of structured events or sequences, as in speech or music”, and that “one is compelled to presume that the power spectrum is an inadequate descriptor, as it would be for speech or music”. We speculate that scale-free brain activity most likely plays an important role in the brain’s operation, and an important goal for future work is to crack the code contained therein. We would not be surprised if new mathematical tools will need to be developed before this goal is within reach. The present study, by showing that across different systems a similar 1/f^β power spectrum can mask unique nested-frequency patterns specific to a system, constitutes an initial, modest step toward this goal.

Linkage to Previous Work

This study follows from a recent line of research emphasizing the physiological significance of broadband LFP, EEG and ECoG activity (Belitski et al., 2008; Bressler et al., 1993; Henrie and Shapley, 2005; Manning et al., 2009; Miller et al., 2009a; Miller et al., 2009b; Montemurro et al., 2008; Onton and Makeig, 2009; Siegel and Konig, 2003; Young et al., 1992). In light of the present results, broadband phenomena could (at least in many cases) be re-interpreted as changes in arrhythmic, scale-free brain activity, which are only sometimes accompanied by brain oscillations (see, e.g., Fig. 4D herein and Henrie and Shapley, 2005; Siegel and Konig, 2003). The recently developed information-theoretic approach may be a promising tool for discerning the functional significance of broadband brain activity (Belitski et al., 2008; Montemurro et al., 2008; Panzeri et al., 2010).

Miller et al. (2009a) described an overall shift of power spectrum in >25 Hz range during a hand motor task, without any significant change in the power-law exponent. The conservation of the power-law exponent in the higher-frequency range is consistent with our finding (Figs. 4 & S4). By contrast, we did not find any overall increase of power in this frequency range. This is likely because we used a widely spaced task paradigm that is much less intensive than the paradigms used in Miller et al. (2009a) and other previous studies showing broadband power increase (e.g., Belitski et al., 2008; Henrie and Shapley, 2005; Manning et al., 2009; Siegel and Konig, 2003; Young et al., 1992). Our paradigm, especially the “cued” condition, is similar to a vigilance task, during which the low-frequency power/variance is reduced (Schroeder and Lakatos, 2009), which accounts for the decrease of power-law exponent that we observed (Figs. 4 & S4). Importantly, this does not mean that the SCP in this frequency range is suppressed. On the contrary, in such tasks, a negative shift of the SCP is expected to develop in activated brain regions, including the “Bereitschaft potential” in the motor cortex and the “contingent negative variation” in the anterior cingulate cortex, corresponding to motor performance and anticipatory attention respectively (Kornhuber and Deecke, 1965; Walter et al., 1964). Similar negative SCP shift can occur in any activated brain region with depolarization of superficial layers (He and Raichle, 2009). This negative SCP shift might be accompanied by a decrease of the fluctuation around the mean, giving rise to the decreased amount of power in the low-frequency range and thus the decreased power-law exponent.

Lastly, the relationship between scale-free brain activity and neuronal avalanches (Petermann et al., 2009) warrants future investigation. In line with our proposal that merely characterizing the power spectrum is insufficient for investigating scale-free dynamics and one has to dig deeper in search for the fine spatiotemporal structures (Fig. 8), these authors stress that a 1/f^β (β ≈ 2) power spectrum is not sufficient for the emergence of neuronal avalanches. Because neuronal avalanches consist of negative peaks of neural field potentials, and the raw waveform of field potentials is the scale-free brain activity giving rise to the 1/f^β power spectrum, a relationship between neuronal avalanches and scale-free brain activity is not only conceivable but indeed seems plausible. Since neuronal avalanches are specific to superficial layers of the cortex (Stewart and Plenz, 2006), a first step might be to investigate whether the 1/f^β power spectrum has any heterogeneity across cortical layers. Surprisingly, we are not aware of any published record on this, and urge such a systematic investigation, which would not only shed light onto the potential relationship between scale-free brain activity and neuronal avalanches, but would also be informative toward the mechanisms underlying scale-free brain activity itself.

The Genesis of Scale-Free Brain Activity

The power-law distribution is a common attractor distribution that can emerge out of many different generative mechanisms (Frank, 2009), and a variety of models have been built to simulate scale-free dynamics. Here we focus on potential neurobiological mechanisms that might underlie scale-free brain field potentials. One obvious candidate is the 1/f^β noise present in neurotransmitter release (Lowen et al., 1997) and neuronal spike trains (Gisiger, 2001; Gruneis et al., 1989; Takahashi et al., 2004; Yamamoto, 1991). However, since these phenomena remain isolated observations, they might not be sufficient to account for the ubiquitous scale-free dynamics in EEG or LFP recordings, and a brain-network perspective will likely also be important.

Notably, functional brain networks extracted from voxel-by-voxel spontaneous fMRI signal correlations have a connectivity distribution that follows power-law scaling with an exponent close to 2, suggesting a scale-free network topology (Eguiluz et al., 2005; van den Heuvel et al., 2008). Furthermore, the geometry of axonal and dendritic trees exhibits self-similarity and scale-invariance properties (Bok, 1959; Freeman, 2007), and simulations show that if axon diameter is dictated by the need to transport materials to the synaptic terminals they would be expected to follow a power-law distribution with an exponent of 2 (Chklovskii and Stepanyants, 2003). While an important future direction in complexity science concerns the bridging of complex networks with complex dynamics (Barabasi, 2009) and the relation between scale-free and small-world networks remains an area of active research (Bullmore and Sporns, 2009), we here tentatively propose that scale-free brain field potentials could be produced by Poisson-like spike trains propagating through scale-free brain networks composed of heterogeneous excitatory neurons, and that the addition of inhibitory neurons introduces a time constant to the network and thus produces periodic brain oscillations added to the scale-free brain activity (another mechanism of producing oscillations is via cellular pacemakers). In line with recent neurophysiological observations (Logothetis et al., 2007), this mechanism does not involve frequency-dependent cortical tissue filtering of field potentials. While our proposal remains to be tested explicitly, it is consistent with recent results from optogenetic manipulations showing that activation of fast-spiking interneurons selectively amplifies gamma oscillations, whereas activation of pyramidal neurons amplifies lower frequencies, the profile of which resembled 1/f^β power spectrum (see Fig. 3d in Cardin et al., 2009). This proposal is also consistent with recent simulations of brain networks (Freeman and Zhai, 2009; Thivierge and Cisek, 2008). In Freeman and Zhai (2009), Poisson-like spike trains were convolved with an impulse response function representing dendritic response and summed, the result of which was scale-free dynamics similar to that recorded by ECoG; the addition of inhibitory neurons introduced narrow-band oscillations that appeared as local peaks rising above the 1/f^β distribution in the power spectrum. Similarly, in Thivierge and Cisek (2008), arrhythmic network spikes were shown to be a direct result of heterogeneity among modeled pyramidal neurons, and the addition of inhibitory pathways introduced higher-frequency activities. This model also showed that global synchronization can emerge out of arrhythmic neuronal activity and that neuronal spikes have long memory – a hallmark of scale-free dynamics.

An influential model of “1/f noise” is the self-organized criticality (SOC) theory (Bak, 1996; Bak et al., 1987), which has been applied extensively to simulate brain networks (Shin and Kim, 2006), neuronal avalanches and “1/f noise” (de Arcangelis et al., 2006; Levina et al., 2007; Lin and Chen, 2005). Although the SOC theory might be more suited to describe neuronal avalanches (Petermann et al., 2009), in our view, convincing evidence suggesting that it is the underlying mechanism giving rise to the 1/f^β activity in raw EEG/LFP activity is still lacking. For example, the SOC theory emphasizes the power-law exponent being close to 1, whereas the power-law exponent of empirical electrophysiological recordings is usually close to 2 (see Fig. 1 herein and Freeman and Zhai, 2009; Miller et al., 2009a; Milstein et al., 2009). Nonetheless, recent extensions of the SOC theory to include a broader range of power-law exponents (De Los Rios and Zhang, 1999) might aid in expanding its explanatory power.

Implications for Psychology and Cognitive Sciences

Scale-free dynamics are widely present in the fluctuations of human cognitive and behavioral performance (Gilden, 2001) and even flight behaviors of Drosophila (Maye et al., 2007). Because the brain evolved through organisms’ struggles in coping with the external world and its main function is to proactively act upon the world through its sense and motor organs for the sake of the organism’s survival, it is hardly surprising that the statistical properties of the brain’s dynamics reflect the statistical properties of the environment and the universally present scaleinvariance thereof; and in turn it is hardly surprising that the cognition and behavior as well as music produced by the human brain follow the same statistical properties. Indeed, several well established psychophysical laws applicable across domains and species, including Weber’s law, can be derived simply from a scale-invariance principle (Buzsaki, 2006; Chater and Brown, 1999). In line with our results showing that the power-law exponent decreased during task performance in activated brain regions, it was found that the slope of the 1/f^β power spectrum in reaction time fluctuations varied parametrically with task difficulty, being steepest for easiest tasks and progressively shallower as task difficulty increased (Ward, 2002). Lastly, given that the power-law exponent of spontaneous ECoG signals is close to 2, it is of interest to note that the power spectra of some cartoons are approximately 1/f¹, whereas those of realistic paintings and photographs are close to 1/f², and impressionist paintings have spectra somewhere in between (Balboa and Grzywacz, 2003; Ward, 2002).

ECoG Data Collection and Power Spectral Calculation

Eight patients undergoing surgical treatment for intractable epilepsy participated in the study. To localize epileptogenic zones, patients underwent a craniotomy for subdural placement of electrode grids and strips followed by 1-2 weeks of continuous video and ECoG monitoring. The placement of electrodes and the duration of monitoring were determined solely by clinical considerations. All patients gave informed consent. ECoG signals were split and sent to both the clinical EEG system and a research EEG system (SynAmp² RT, Neuroscan, DC-coupled recording). All data in the present study were from the research amplifier. For details on patient and data collection information see Supplementary Experimental Procedures and Supplementary Table 1.

fMRI Data Collection and Analyses

Resting-state fMRI data (with visual fixation) were collected from 17 normal healthy volunteers on a 3T Siemens Allegra MR scanner. Each subject completed 4 fMRI runs, each run lasting 7 min. Standard fMRI data pre-processing procedure was used. The normalized power spectra were derived from the lagged auto-correlation functions using a Bartlett smoothing procedure. For details see Supplementary Experimental Procedures.

Earth Seismic Wave and Stock Market Data Collection

Continuous natural seismic wave data (sampling rate 100 Hz) were collected by University of Nevada, Reno seismic network. The historical daily prices of Dow-Jones Industrial Average index from Oct 1, 1929 to Apr 9, 2009 were downloaded from Yahoo Finance website.

Nested-Frequency Analysis

To assess cross-frequency phase-amplitude coupling, an inverted entropy measure was used to yield a modulation index (MI), following Tort et al. (2008). The MI was then compared with shuffled data to obtain a Z-score measurement of cross-frequency coupling strength (MI Z-score). For details on the computation of the MI and the shuffling procedure see Supplementary Experimental Procedures.