Abstract
A growing body of evidence suggests that the neuronal dynamics are poised at criticality. Neuronal avalanches and long-range temporal correlations (LRTCs) are hallmarks of such critical dynamics in neuronal activity and occur at fast (subsecond) and slow (seconds to hours) timescales, respectively. The critical dynamics at different timescales can be characterized by their power-law scaling exponents. However, insight into the avalanche dynamics and LRTCs in the human brain has been largely obtained with sensor-level MEG and EEG recordings, which yield only limited anatomical insight and results confounded by signal mixing. We investigated here the relationship between the human neuronal dynamics at fast and slow timescales using both source-reconstructed MEG and intracranial stereotactical electroencephalography (SEEG). Both MEG and SEEG revealed avalanche dynamics that were characterized parameter-dependently by power-law or truncated-power-law size distributions. Both methods also revealed robust LRTCs throughout the neocortex with distinct scaling exponents in different functional brain systems and frequency bands. The exponents of power-law regimen neuronal avalanches and LRTCs were strongly correlated across subjects. Qualitatively similar power-law correlations were also observed in surrogate data without spatial correlations but with scaling exponents distinct from those of original data. Furthermore, we found that LRTCs in the autonomous nervous system, as indexed by heart-rate variability, were correlated in a complex manner with cortical neuronal avalanches and LRTCs in MEG but not SEEG. These scalp and intracranial data hence show that power-law scaling behavior is a pervasive but neuroanatomically inhomogeneous property of neuronal dynamics in central and autonomous nervous systems.
Introduction
The human brain appears to operate akin to a dynamical system near a critical point (Chialvo, 2010), which could be central for understanding neuronal processing in vivo (Shew et al., 2009, 2011) and its implications for behavior (Palva et al., 2013) and brain diseases. Critical dynamics is an emergent characteristic also in several models of brain activity (Rubinov et al., 2011; Poil et al., 2012). Nevertheless, the predictions from theories and models of criticality are challenging to test with empirical data from biological systems (Mora and Bialek, 2011). For criticality in brain dynamics, two principal series of experimental observations constitute the most compelling evidence.
First, in the timescales of seconds to tens of minutes, power-law long-range temporal correlations (LRTCs) characterize the amplitude envelopes of neuronal oscillations in human MEG and EEG (Linkenkaer-Hansen et al., 2001) as well as in intracranial recordings (Monto et al., 2007). The LRTC scaling exponent, β, is an individual characteristic that is modulated by stimuli and tasks and is different among frequency bands (Linkenkaer-Hansen et al., 2004; Palva et al., 2013). Furthermore, fMRI has shown that the resting-state LRTC exponents of BOLD signals are also different among brain systems (He, 2011), but this has not been systematically assessed with electrophysiological recordings so far.
Second, in millisecond timescales, neuronal activity cascades, or “neuronal avalanches,” exhibit power-law size distributions in in vitro (Beggs and Plenz, 2003; Pasquale et al., 2008), in vivo (Petermann et al., 2009; Hahn et al., 2010), and human invasive (Solovey et al., 2012; Priesemann et al., 2013) and noninvasive (Palva et al., 2013; Shriki et al., 2013) electrophysiological recordings. Neuronal avalanches in organotypic cultures (Beggs and Plenz, 2003) and human MEG (Shriki et al., 2013) appear to exhibit balanced branching at the size scaling exponent α of 3/2, which is line with a theoretical prediction for a critical branching process (Zapperi et al., 1995).
Several studies show that neuronal avalanches may coexist with LRTCs both in vitro (Beggs and Plenz, 2004) and in vivo (Benayoun et al., 2010; Hahn et al., 2010; Palva et al., 2013) and in models (Poil et al., 2012). LRTCs and avalanches may hence arise from shared underlying mechanisms or coemerge through interactions with a third system, such as the autonomous nervous system (ANS) that also exhibits scale-free activity (Ivanov et al., 1999).
We use here both intracranial stereo-electroencephalography (SEEG) and source-reconstructed MEG to obtain, at mesoscopic and macroscopic, respectively, levels, a comprehensive view onto the electrophysiological hallmarks of criticality in human resting-state brain activity. We set out here to do the following: (1) rigorously assess the scaling and branching behavior of neuronal avalanches; (2) characterize the cortical topology of LRTCs and its relation to brain systems in fMRI; and (3) test the hypothesized relationship between LRTCs and avalanches. Finally, we also addressed the relationship of criticality in the central nervous system with the scaling of heart-rate variability (HRV) that is controlled by the ANS.
Materials and Methods
Data acquisition.
We analyzed the invasive SEEG recordings from a cohort of 22 epileptic patients (7 females, age 16–21 years) and noninvasive MEG data recorded from 14 healthy subjects (7 females, age 18–27 years). Resting-state SEEG was collected for 10 min with eyes closed and without external disturbance using a 192 channel SEEG amplifier system (NIHON-KOHDEN NEUROFAX-110) at a sampling rate of 1 kHz. We acquired monopolar local field potentials (LFPs) from brain tissue with platinum-iridium, multilead electrodes. The number of electrode contacts along each penetrating shaft varied from 8 to 15. These contacts were 2 mm long, 0.8 mm thick, and had an intercontact distance of 1.5 mm (DIXI Medical). The anatomical positions and amounts of electrodes varied according to surgical requirements (Cardinale et al., 2013). On average, each subject had 14 ± 1.9 (mean ± SD) shafts (range 17–10) with a total of 152 ± 20 electrode contacts (range 193–147, left hemisphere: 37 ± 49, right hemisphere: 115 ± 51 contacts). MEG data were collected in resting-state condition for 10 min with the subjects looking at a fixation point on the monitor screen. We recorded 306 channel (204 planar gradiometers and 102 magnetometers) MEG (Elekta Neuromag) at a sampling rate of 600 Hz.
SEEG preprocessing.
The locations of each SEEG electrode contact were extracted with submillimeter accuracy from postimplant cone-beam CT scans by means of an automatic algorithm for SEEG implant identification (Arnulfo et al., 2015a) and subsequently coregistered to FreeSurfer (http://surfer.nmr.mgh.harvard.edu/) geometrical space. We used a “closest-white-matter” referencing scheme for SEEG where electrode contacts in gray matter were referenced to the closest contacts in underlying white matter (Arnulfo et al., 2015b), which largely eliminates volume-conducted signals but preserves the LFP phase and polarity undistorted because gray matter signals are never referenced to other gray matter signals. Only signals from gray-matter contacts were analyzed in this study.
MEG preprocessing and source reconstruction.
Preprocessing of raw MEG time series consisted of the following steps. First, the temporal extension of signal space separation method (Taulu et al., 2005) was used to remove extracranial noise from the raw MEG recordings. Next, independent component analysis (Bell and Sejnowski, 1995) was used to identify and exclude components associated with eyes movements/blinks and cardiac artifacts.
For cortical surface reconstructions, T1-weighted anatomical magnetic resonance images with a 1.5-T MRI scanner (Siemens) were recorded. FreeSurfer software was used for automatic volumetric segmentation of the MRI data, surface reconstruction, flattening, and cortical parcellation. MNE software (www.martinos.org/mne/) was used to create three-layer boundary element conductivity models and cortically constrained source models with fixed-orientation dipoles, and for computing the forward and inverse operators (Hämäläinen and Ilmoniemi, 1994; Lin et al., 2006). MEG time series were filtered using filter with a finite impulse response (pass-band 1–40 Hz). After filtering, MEG sensor data were inverse transformed and then collapsed into time series of 219 cortical parcels derived from individual MRI scans using FreeSurfer and CMP software (Daducci et al., 2012) using collapsing operators maximizing individual reconstruction accuracy (Korhonen et al., 2014).
Analysis of neuronal avalanches.
We identified neuronal avalanches from source-reconstructed MEG parcel time series and from white-matter referenced SEEG gray-matter LFP time series and quantified their statistical properties in the following manner. Broadband-filtered time series (1–40 Hz) were normalized by subtracting the mean and dividing by SD. The time series were transformed into binary point process by detecting peaks above the threshold T that ranged from 1.5 to 5.25 with step 0.25 and setting the samples corresponding to peak latencies to 1's. These binary sequences (or sequences of events) were then converted into avalanche time series by summing the events across the electrodes in different time bins (Δt, ranges from 4 to 80 ms with step 4 ms) (Plenz, 2012). A neuronal avalanche is defined as cluster of events in successive time bins where the beginning and end of the avalanche are defined by single time bins with no events (see Fig. 1A). The avalanche size is the total number of events. It has been shown earlier that the avalanche size distribution is typically a power-law (Beggs and Plenz, 2003; Plenz, 2012; Palva et al., 2013; Shriki et al., 2013). Importantly, the distribution reveals power-law form only for a certain parameters (threshold and time bin width) used in avalanches detection. To comprehensively map the parameter space, we assessed the avalanche sizes for multiple combinations thresholds and time bin widths and fit the data using the truncated-power-law model, which has separate power-law and exponential components. The truncated-power-law model was statistically tested (see below) against power-law and exponential to define the range of parameters (T and Δt) or regimen where the one model outperforms another.
Probability functions of truncated-power-law, power-law, and exponential models can be expressed as follows:
where p(s) denotes the probability of observing avalanche of size s, Cp, Ctp, and Ce are normalization constants, α and λ are parameters of the power-law and exponential terms, and smin and smax are the minimum and maximum avalanche sizes, respectively.
There are several studies that emphasize the importance of estimating minimum and maximum avalanche sizes to obtain correct values of scaling exponents (Clauset et al., 2009; Klaus et al., 2011; Dehghani et al., 2012; Langlois et al., 2014). The minimum size is highly affected by noise level, which depends on the thresholds and time bin widths used in avalanche detection. Instead of estimating smin, we used fixed value (smin = 1) and assessed the scaling exponents for multiple thresholds and time bin widths; thus, the exponent α can be considered as a function of two parameters (T and Δt), α(T,Δt). The maximum avalanche size smax was set to the total number of cortical parcels (MEG) or electrode contacts (SEEG) included in the analysis (Priesemann et al., 2013; Shriki et al., 2013).
To estimate the scaling exponents of avalanche size distribution, maximum likelihood method was applied (Clauset et al., 2009; Klaus et al., 2011). The results of maximum likelihood approach were corroborated with a freely available Python toolbox (Alstott et al., 2014) for scaling exponent estimation (Clauset et al., 2009).
Definition of the regimens.
We used statistical testing to define three regimens that correspond to a certain combinations of T and Δt. Power-law regimen assumes no difference between truncated-power-law and power-law models, as well as exponential regimen that assumes no difference between truncated-power-law and exponential models. Truncated-power-law regimen corresponds to the situation where power-law or exponential models are different from truncated-power-law.
Statistical comparison between the models has been done by computing log-likelihood ratio (LLR) as it described in (Clauset et al., 2009; Klaus et al., 2011).
The LLR between the two distributions is defined as follows:
where a and b are estimated parameters of the first and second model. If LLR(x) is significantly larger than 0, then the first model is assumed to be better fit data x, and vice versa. The p value for the LLR test is defined as follows:
where:
Difference between two models was considered to be statistically significant if p < 0.05.
Likelihoods for truncated-power-law, power-law, and exponential models can be expressed as follows:
Analysis of long-range temporal correlations.
We used detrended fluctuation analysis (DFA) to assess the scaling exponents of LRTCs (Peng et al., 1995; Linkenkaer-Hansen et al., 2001; Palva et al., 2013). DFA was applied to the amplitude envelopes of neuronal time series that are filtered using Morlet's wavelets with the logarithmically spaced central frequencies (from 3 to 40 Hz). The same analysis was used to assess LRTCs of heartbeat intervals of MEG subjects. The analysis can be represented as a two stage procedure. In the first stage, time series X is normalized to 0 mean and integrated over the samples (y(n) = (X(n) − 〈X〉) + Y(n − 1)), then segmented into time windows of various sizes τ. At the second stage, each segment of integrated data is locally fitted to a linear function and the mean-squared residual F(τ) is computed,
where N is the number of samples in segment τ.
The scaling exponent β is defined as the slope of linear regression of the function F(τ) plotted in log-log coordinates. The exponents β(f) were computed for each wavelet central frequency and cortical parcel (MEG) or electrode contact (SEEG), in case of neuronal recordings, whereas βHRV was evaluated for the HRV time series extracted from MEG data using independent component analysis.
Relationship between scaling exponents α and β.
We assessed the relationship between the pairs of scaling exponents of avalanche size distribution (α) and those of LRTCs (β) (associated with each subject) with Pearson correlation coefficient. The significance of the correlation coefficient was assessed with t-statistics, t =
Partial correlation analysis has been applied to control the potential confounders. The partial correlation was computed between α and β with controlling βHRV, and α and βHRV with controlling β.
Shuffling methods.
To assess statistical significance of the results without making an assumption on the theoretical distribution of data, the scaling exponents of LRTCs and neuronal avalanches were computed for phase-shuffled data (Linkenkaer-Hansen et al., 2001; Shriki et al., 2013). Phase-shuffling disrupts temporal as well as spatial correlations in multichannel time series while preserving the power spectrum. Phase-shuffling is often used in hypothesis testing for avalanche size distribution (Gireesh and Plenz, 2008; Shriki et al., 2013). Here the problem is that, although spatial interactions, such as neuronal activity cascades (Plenz, 2012), would appear to be a likely explanation for the power-law statistics of neuronal avalanches in empirical data, similar scaling could arise artificially through the local temporal autocorrelation structures of thresholded long-memory process (Touboul and Destexhe, 2010). Whether the power-law size scaling of avalanches is attributable to global spatial interactions or local temporal correlations cannot be deduced with phase-shuffled surrogate data. For this reason, scaling exponents of avalanche size distribution were assessed also with time-shuffled surrogate data. The time-shuffled time series were derived from the original ones so that the signal of each channel was circularly shifted by a random number of samples. This approach breaks spatial correlations but essentially fully preserves the temporal autocorrelation structure.
Estimation of the branching parameter.
The branching parameter (σ) was estimated by computing the ratio of the number of events in the second time bin to the number of events in the first time bin in neuronal avalanche (Beggs and Plenz, 2003).
Where ni(1) and ni(2) is the number of events in first and second time bin of ith avalanche, respectively; N is the total number of avalanches.
Functional parcellation.
To assess the putative anatomical differences in neuronal scaling properties with a readily accessible neuroanatomical basis, we collapsed the MEG data and the large amount of SEEG recording sites into functional modules that are consistent across subjects (Damoiseaux et al., 2006; Brookes et al., 2011). The scaling exponents of LRTCs, β(f), and correlation coefficients, r(β(f), α(T, Δt)) and r(β(f), βHRV), were hence mapped into fMRI-based functional systems (Yeo et al., 2011). Each cortical parcel (MEG) or electrode contact (SEEG) was associated with one of the seven functional systems, including visual (VI), somato-motor (SM), dorsal attention (DA), ventral attention (VA), limbic (LI), frontoparietal (FP), and default (DM). The mapping was done for each subject by comparing 3D coordinates of the contacts or source points within parcel with the coordinates of functional systems and assigning the parcel or contact with closest (smallest Euclidian distance) system. Consensus clustering was applied to unify the functional parcellation between the subjects.
Results
Avalanche size distributions reveal scale-free dynamics of neuronal fluctuations
Cascades of propagating neuronal activity, neuronal avalanches, were readily observable in broadband source-reconstructed MEG (Fig. 1A) as well as in gray-matter SEEG recordings of spontaneous activity. The avalanche is defined to comprise consecutive events that for field potential data are peaks with amplitude greater than T in neuronal fluctuations, summed across cortical parcels (MEG) or electrodes (SEEG) within a time bin of a certain width (Δt) (Fig. 1A). The avalanche size is defined as the number of events on all sensors in the avalanche (Shriki et al., 2013) (Fig. 1A). We assessed the scaling exponents of avalanche size distributions, α, by using the maximum likelihood approach to compare power-law, truncated -power-law, and exponential models.
Large-magnitude neuronal fluctuations constitute neuronal avalanches in MEG and SEEG. A, Example of spatiotemporal pattern of neuronal avalanches at fast timescales. Avalanche time series (black bars) are derived from multichannel recordings (color lines), where the peaks (black circles) with suprathreshold amplitude (T, green line) are transformed to binary 1's and summed across channels within the time bins (Δt, vertical lines). Neuronal avalanches are described by their size (s; i.e., number of peaks the cascade contains). B, Example of a parametric map of scaling exponents of avalanche size distribution obtained for multiple thresholds (T) and time bin widths (Δt). C–E, Example of distributions of avalanche sizes detected at high (C), moderate (D), and low (E) thresholds obeying power-law, truncated-power-law, and exponential form, respectively. Colored circles represent the parameters (T and Δt) for which neuronal avalanches are detected. Gray lines/circles represent distribution for phase-shuffled time series.
The estimates of α as well as the best-fitting model per se are dependent on T and Δt, which may lead to controversial results if a priori fixed parameter values are used in data analysis. To quantify α comprehensively, we explored a wide range of T and Δt (Fig. 1B). In both MEG and SEEG, we found three salient regimens in these data: at high thresholds, the size distribution was well fit by a power-law (Fig. 1C), whereas at intermediate and low thresholds, the distributions were best fit by a truncated power-law or an exponential, respectively (Fig. 1D,E). The fitting range was restricted by the total number of cortical parcels (MEG) or electrode contacts (SEEG), as has been suggested previously (Priesemann et al., 2013; Shriki et al., 2013). The avalanche-like events detected with an identical procedure for phase-shuffled time series were not well fit by the power-law model with any parameters.
Relationship of the scaling exponent and branching parameter of neuronal avalanches
Scaling exponents assessed for multiple thresholds and time bin widths were averaged across subjects. The regimens (ranges of the parameters T and Δt) indicating significant difference between truncated-power-law, power-law, and exponential models were defined at the group level assuming that the individual models belong to the same class and significantly different (log-likelihood test) from alternative models for given T and Δt (see Materials and Methods).
We found a systematic increase in the exponent values with increasing threshold and decreasing time bin width for original but also time-shuffled MEG and SEEG time series (Fig. 2A). The scaling exponent α can be expressed as a function of T and Δt (α(T,Δt) = Δt−k/T) where k is a normalizing parameter. This observation was in line with previous findings where α is considered as a function of Δt for relatively high thresholds (Beggs and Plenz, 2003; Petermann et al., 2009). The results suggested that permanent value of α (e.g., α = 3/2) requires simultaneous changing of the threshold and time bin width, which is in line with the limited range of Δt where fixed T yields unchanged α (Shriki et al., 2013).
Scaling exponents of avalanche size distribution vary as a function of threshold and time bin width and bounded by power-law and exponential regimens. A, Scaling exponents computed for multiple thresholds and time bins (averaged across subjects) show a clear dependency on threshold and time bin widths for original (left) and time-shuffled (right) data. Exponents that are close to those of a critical branching process (α(T,Δt) = 3/2) are approximated by a power-law function of threshold and time bin width (black line with yellow circles). Statistical estimates of power-law (gray lines) and exponential (yellow lines) regimens reveal different ranges of thresholds and time bins for MEG (top) and SEEG (bottom) time series while showing a similar trend. B, Dynamics of the branching parameters is similar to that of scaling exponents across multiple thresholds and time bin widths for original (left) and time-shuffled (right) data. The branching parameters that are close to those of a critical branching process (σ = 1.0) are approximated by a power-law function of T and Δt (black line with cyan circles) for MEG (top) and SEEG (bottom) time series. Yellow circles represent the scaling exponents close to 3/2. Gray lines indicate the power-law regimen. Yellow lines indicate exponential regimens. Scaling exponents for a certain combinations of T and Δt are shown in dark gray because of insufficient amount of data for statistical analysis or the values are out of range.
The exponents computed with the same T and Δt were larger for SEEG than for MEG time series (p < 0.0005, t test). This observation may be explained by linear mixing in MEG and the exponents being negatively correlated with linear mixing (Shriki et al., 2013) or by the putative anatomical subsampling of SEEG (Priesemann et al., 2009). In both SEEG and source-reconstructed MEG, the scaling exponents at the commonly used parameters T = 3 SD and Δt = 4 ms were greater than those reported earlier for sensor space MEG time series (Fig. 2A) (Shriki et al., 2013). Avalanches in phase-shuffled surrogate data did not exhibit power-law size scaling, as observed previously (Shriki et al., 2013). However, avalanches in time-shuffled surrogate data, where spatial structures were abolished but the temporal autocorrelation structure retained, did exhibit power-law size distributions. Power-law scaling in time-shuffled data was, nevertheless, associated with greater scaling exponents (p < 0.05, t test, FDR corrected) in the power-law regimen than the original data in >85% of all analyzed pairs of T and Δt. This suggests that power-law size scaling of neuronal avalanches does not arise epiphenomenally from the spurious superpositions of local long-memory random processes (Bédard et al., 2006; Touboul and Destexhe, 2010).
To test the hypothesis that the neuronal dynamics in MEG and SEEG time series was akin to a critical branching process (Plenz, 2012), we computed the branching parameter (σ) for multiple T and Δt in the same manner as the scaling exponents (Fig. 2B). The results showed that the values of σ observed within the power-law regimen of the T-Δt plane were smaller than those expected for a critical branching process (σ = 1). Nevertheless, the curves α = 3/2 and σ = 1 were very close to each other and contained within the truncated-power-law regimen (Fig. 2B). Branching parameter for time-shuffled data revealed a similar trend without a significant difference (p > 0.11, MEG; p > 0.10, SEEG; ANOVA) between the parameters that yielded α = 3/2 and σ = 1 (Fig. 2B). These surrogate data hence show that the temporal autocorrelation structure per se in MEG and SEEG recordings of large-scale neuronal network activity, in the absence of true spatial correlations, leads to similar truncated-power-law scaling of avalanche sizes with seemingly balanced branching. This result questions the validity of using the theory of branching process in describing neuronal dynamics in modular large-scale systems, which has also been discussed previously (Taylor et al., 2013; Hartley et al., 2014).
Scaling exponents of LRTCs reveal signatures of functional and frequency specificity
We then quantified the neuronal LRTCs with scaling exponents β that were estimated for the amplitude envelopes of narrowband oscillations (Fig. 3A) across a wide range of frequencies. Colocalizing the electrophysiological data with a functional-connectivity-fMRI (fc-fMRI) based parcellation of intrinsic connectivity networks, we asked whether LRTCs were distinct among functional brain systems (Fig. 3B). The exponents of LRTCs were assessed using DFA that quantifies the scaling of the signal autocorrelation structure of time series through computing variance of linearly detrended signal at different timescales. The DFA plots were log-log linear across a wide time range (from 5 to 500 s), which implies monofractality and thereby that a single exponent β is sufficient to describe dynamics of the process (Fig. 3C). The scaling exponents were largely above β = 0.5 that would correspond to an uncorrelated Gaussian process, and varied along the frequency of underlying narrowband oscillations, revealing a prominent peak at 10 Hz in MEG and less pronounced at 8 Hz in SEEG (Fig. 3D). The exponents associated with five frequency bands (δ-band (3 Hz), θ-band (6 Hz), α-band (10 Hz), β-band (16 Hz), and γ-band (30 Hz)) were mapped into seven functional systems, including VI, SM, DA, VA, LI, FP, and DM (Fig. 3E). The power-law exponents β(f) ranged from 0.61 (LI, in δ-band) to 0.72/0.70 (VI/DA, in α-band) in MEG and from 0.63 (LI, in β-band) to 0.70 (DA, in α-band) in SEEG. The results showed a good agreement between MEG and SEEG (Pearson correlation coefficient r = 0.39, p < 0.02; and r = 0.49, p < 0.005), with and without VI system, respectively) with the VI system as an outlier possibly because of sparse SEEG sampling there. We then assessed the similarity of β(f) between the systems separately for MEG and SEEG.
Scaling exponents of LRTCs are unique characteristics of cortical oscillation amplitude dynamics. A, Example of narrowband filtered source-reconstructed time series (gray line) of a single cortical parcel and its envelope (brown line) for which LRTCs are computed. B, Seven fMRI-based functional systems are shown on an inflated brain surface: visual (blue), somatomotor (cyan), dorsal attention (green), ventral attention (teal), limbic (brown), frontoparietal (yellow), and default (red). C, DFA reveals mono-fractal behavior of LRTCs (black circles/line) in MEG (top) and SEEG (bottom) time series of representative subjects (f = 10 Hz). The exponents for phase-shuffled time series (βref) are close to those of random Gaussian process (gray line). D, Scaling exponents (pooled across sensors/subjects) vary with the frequency of underlying narrowband oscillations. The mean values of scaling exponents (averaged across sensors/subjects) as a function of frequency show maximal values at frequencies ∼10 Hz and 8 Hz (black line) for MEG and SEEG data, respectively. Phase-shuffled data are characterized by smaller exponents and frequency-unspecific dynamics (gray lines). E, The exponents of LRTCs are significantly different (p < 0.001, FDR corrected, t test) between some of the functional systems at certain frequencies. Bar colors match functional system colors as shown in panel B.
We found that, in MEG, the exponents in θ-, α-, and β-frequency bands were different in VI, SM, and DA compared with LI, FP, and DM (Fig. 3E; p < 0.001, unequal variance t test). In the γ-frequency band, differences were observed between VA and DM, and SM and DA. Importantly, the lowest values of β(f) were observed in LI in MEG and SEEG in θ-, α-, and β-bands.
These data indicate that, in the light of both mesoscopic and macroscopic electrophysiological measurements, spontaneous resting-state brain activity exhibits significant LRTCs throughout the human neocortex with distinct scaling properties in different brain systems and frequency bands.
Scaling exponents of avalanche size distributions are correlated with those of LRTCs
We characterized the relationship between the scaling exponents of avalanche size distributions and LRTCs. We first averaged the exponents β(f) across cortical parcels (MEG) or electrode contacts (SEEG) for each subject and estimated the correlation between β(f) and α(T,Δt) for parameters of T and Δt picked from the power-law and truncated-power-law regimens. Surprisingly, we observed strong negative correlations (r = −0.78, p < 0.005, MEG; r = −0.44, p < 0.05, SEEG) for parameters (f = 3 Hz, T = 3.25, Δt = 8, MEG; and Δt = 16, SEEG) that belong to power-law regimen and positive correlations (r = 0.80, p < 0.005, MEG; r = 0.46, p < 0.05, SEEG) for parameters for the truncated-power-law regimen (f = 3 Hz, T = 2.00, Δt = 8, MEG; and Δt = 16, SEEG) (Fig. 4A). Despite the fact that correlation coefficients were higher for MEG than SEEG, the results were strikingly consistent between these modalities. To assess this phenomenon comprehensively, the correlations were computed for all T and Δt in the avalanche analysis, and for each frequency of LRTCs (f), and then significant values were averaged across frequencies (Fig. 4B). The results indicated negative correlations that broadly represented within power-law regimen, whereas positive correlations were observed within a narrow range of parameters (T and Δt) that precisely coincide with α(T,Δt) = 3/2 in truncated-power-law regimen. The values of negative and positive correlations were depended on frequency of β(f) and varied independently from each other in power-law and trunctated-power-law regimens.
Scaling exponents (α) of avalanche size distributions are correlated with LRTC exponents (β). A, Scatter plots of α(T,Δt) and β(f) for all subjects for avalanches belonging to the power-law regimen (f = 3 Hz, T = 3.25 SD, Δt = 8 ms [MEG] and Δt = 16 ms [SEEG]) indicate strong negative correlations (blue dots), whereas the scatter plots for avalanches belonging to the truncated-power-law regimen (f = 3 Hz, T = 2.00 SD, Δt = 8 ms [MEG] and Δt = 16 ms [SEEG]) reveal positive correlations (red dots), for MEG (top) and SEEG (bottom) time series. **p < 0.005. *p < 0.05. B, Correlation maps between scaling exponents (α(T,Δt) and β(f)) estimated for multiple thresholds and time bins (significant values are averaged across frequencies) reveal robust patterns of negative values that belong to a power-law regimen (gray lines) and positive correlations that coincide with α(T,Δt) = 3/2 (black line with yellow circles). C, D, Example of avalanche size distributions of two representative subjects with weak (orange) and strong (black) LRTCs for truncated-power-law (C) and power-law (D) regimens, in MEG (top) and SEEG (bottom). Color codes correspond to those in A.
To understand the possible mechanism of interaction between α(T,Δt) and β(f), we assessed the distribution of avalanche size for subjects with weak versus strong LRTCs (Fig. 4C,D). Subjects with weaker LRTCs had smaller exponents α at low thresholds and larger at high thresholds compared with subjects with stronger LRTCs, which might explain the presence of negative and positive correlations between α(T,Δt) and β(f).
Surprisingly, a qualitatively similar relationship was observed for time-shuffled data (Fig. 5), albeit with weaker positive correlations than those of original. These positive and negative correlations were largely observed within the truncated-power-law regimen and coincided with α(T,Δt) = 3/2.
Scaling exponents of avalanche size distributions are correlated with the exponents of LRTCs for time-shuffled data. A, Scatter plots of α(T,Δt) and β(f) across subjects for avalanches belonging to a power-law regimen (f = 3 Hz, T = 3.25 SD, Δt = 8 ms [MEG] and Δt = 16 ms [SEEG]; blue dots) indicate strong negative correlations, whereas the scatter plots for avalanches from truncated-power-law regimen (f = 3 Hz, T = 2.00 SD, Δt = 8 ms [MEG] and Δt = 16 ms [SEEG]; red dots) reveal positive correlations, for MEG (top) and SEEG (bottom) time series. **p < 0.005. B, Correlation maps between scaling exponents (α(T,Δt) and β(f)) estimated for multiple thresholds and time bins (significant values are averaged across frequencies) reveal negative and positive values that largely belong to a truncated-power-law regimen and are close to α(T,Δt) = 3/2 (black line with yellow circles). Gray lines indicate power-law regimen. Yellow lines indicate exponential regimen.
These observations hence comprise the first cortex-wide view into the relationship of avalanches and LRTCs. The finding that their correlation was observed in time-shuffled surrogate data as well, albeit weakly, suggests that, in addition to interareal neuronal interactions, also local LRTCs play a significant role in shaping the avalanche dynamics.
Functional localization of correlates between α and β
The localization of the correlates in terms of functional systems and frequencies of LRTCs was done for earlier defined parameters (T and Δt) associated with negative and positive correlations in MEG data. SEEG data have been excluded from the analysis because of spatial sampling bias. The correlation coefficients were averaged within each functional system and frequency band (Fig. 6). The negative and positive correlates revealed distinct dynamics over the frequencies. In all the frequency bands except θ-band, the negative correlations were significant (p < 0.05), whereas positive correlations were significant only at δ- and γ-frequencies (p < 0.05). The negative and positive correlations showed similar relationships between the systems at δ-frequency band where DA reveals the highest value compared with most systems (VI, LI, FP, and DM). Surprisingly, the negative correlations in α-band were lowest in the VI system compared with others (SM, VA, and LI) (Fig. 6A). DA showed dramatic drop in β-band and became smallest compared with VI, VA, FP, and DM. No differences between the systems were observed in γ-band. The LI revealed highest values compared with all the systems in γ- and partially β-bands in case of positive correlates (Fig. 6B). The results suggested that the spatio-spectral patterns of negative and positive correlates are highly overlapping (r = 0.35, p < 0.03). The spatio-spectral patterns of correlations in time-shuffled data (Fig. 6C) were similar to the negative (r = 0.54, p < 0.0007) and positive (r = 0.76, p < 0.0001) patterns of original time series, although no similarities were observed for phase-shuffled data (Fig. 6D). These findings further corroborate the correlational relationship of LRTCs and neuronal avalanches and raise the intriguing question of whether another long-memory process would through c-modulation determine the dynamics of both neuronal LRTCs and avalanches.
Spatio-spectral patterns of correlates (between α(T,Δt) and β(f)) reveal both system- and frequency-band-specific dynamics. A, Negative correlations observed in the power-law regimen (T = 3.25, Δt = 8 ms) are significant for most frequency bands, except θ, and show differences between systems (*p < 0.01, t test). Dashed line indicates level of significance (p < 0.05, t test for correlation coefficient). B, Positive correlations that belong to truncate-power-law regimen (T = 2.00, Δt = 16) are significant at δ and γ bands, with dominance of limbic (LI) system at high frequencies. C, Negative correlations observed in power-law regimen (T = 3.25, Δt = 8 ms) for time-shuffled data are significant at θ and γ frequencies. D, Phase-shuffling data reveal no significant correlations between the scaling exponents. The functional systems are shown in blue (VI), cyan (SM), green (DA), teal (VA), brown (LI), yellow (FP), and red (DM).
Relationships among LRTCs of neuronal and HRV and dynamics of neuronal avalanches
Our previous study revealed that the scaling exponents of LRTCs of neuronal amplitude fluctuations and of HRV time series were correlated (Palva et al., 2013), which suggests that the critical dynamics in the central and ANSs could be related. To expand this observation and assess its role in the LRTC-avalanche relationship examined here, we first quantified the LRTCs in HRV and examined their relationship with neuronal avalanches. Persistent LRTCs were observed in HRV time series (Fig. 7A), with DFA revealing scale-invariant mono-fractal dynamics (Fig. 7B). The scaling exponents βHRV were correlated with the exponents of avalanche size distribution α(T,Δt) detected in MEG in a manner similar to their correlation with β(f), which supports the notion on a strong interconnection between central and ANS activities in a wide range of timescales (Fig. 7C). Moreover, although the negative correlations between βHRV and α(T,Δt) in MEG were observed for higher values of T and Δt compared with β(f), the positive correlations were precisely matched with α(T,Δt) = 3/2 and σ(T,Δt) = 1 in the truncated-power-law regimen, which suggests stability of the dynamics in this regimen. However, the correlations between βHRV and α(T,Δt) were not significant in SEEG data (Fig. 7D).
The LRTC scaling exponents of HRV are correlated with the size-scaling exponents of neuronal avalanches in MEG but not in SEEG. A, Example of a heartbeat intervals time series for a representative MEG subject. B, DFA reveals scale-free dynamics of HRV time series. C, D, Scaling exponents of HRV (βHRV) are correlated with the exponents of avalanche size distribution α(T,Δt) in MEG (C) but not in SEEG (D).
We then assessed the relationship between βHRV and neuronal LRTCs. Joint distribution of the exponents, βHRV and β(f) (f = 3 Hz, averaged across cortical parcels [MEG] or electrode contacts [SEEG]), plotted for all the MEG and SEEG subjects, indicated a strong linear relationship with a correlation coefficient of r = 0.59 (p < 0.05) in MEG (Fig. 8A). The correlations between βHRV and β(f) were frequency-specific changes and different among functional systems (Fig. 8B). The correlations were significant (p < 0.05) in δ- and γ-frequency bands in MEG. LI and VA systems showed greater correlation values than VI, SM, and DM in δ-band. In γ-band, LI showed the highest value as well, which was different from DA, FP, and SM, and also VI revealed difference with SM. Curiously, the relationship between the systems in γ-band was very similar to that of positive correlations between α(T,Δt) and β(f). Assessing this quantitatively, we found that the spatio-spectral correlation patterns of βHRV and β(f) were similar to the correlation patterns of α(T,Δt) and β(f) for original (r = 0.66, p > 0.0001 and r = 0.72, p < 0.0001; negative and positive correlates, respectively) as well as time-shuffled (r = 0.93, p > 0.0001; negative correlates) data.
LRTCs of HRV and neuronal fluctuations are correlated in MEG but not SEEG. A, Scatter plot of LRTC scaling exponents of HRV (βHRV) and neuronal fluctuations (β(f), f = 3 Hz) for MEG (top) and SEEG (bottom) time series. B, Correlations between βHRV and β(f) reveal frequency-specific changes and differences between functional brain systems (*p < 0.01, t test) in MEG (top) but not SEEG (bottom). Dashed line indicates the significant threshold for correlation coefficients (p < 0.05, t test for correlation coefficient). Functional system colors are the same as in Figures 3 and 6.
In SEEG, however, similarly to the absence of correlations between βHRV and α(T,Δt), there were no significant correlations between βHRV and β(f) in any functional system or frequency band (Fig. 8B). Correlations for some systems seemed relatively high (such as VI at 6 Hz, r = −0.44; LI at 6 Hz, r = −0.45), but only half of the subjects (11 of 22) had electrode contacts in these systems; therefore, the statistical power may be inadequate. Nevertheless, this result strongly suggests that neuronal activity at the scale of millimeter-range LFPs is essentially independent of ANS dynamics that, in turn, is correlated only with the centimeter-scale local synchronization measured by MEG.
Factors determining the correlation between scaling exponents
We applied partial correlation analysis to control the influence of a third variable among the interactions of α(T,Δt), β(f), and βHRV. Intriguingly, we found that the scaling exponents of neuronal avalanches were independently correlated with either cortical (Fig. 9A) or autonomous (Fig. 9B) nervous system LRTCs. Hence, neither central nor ANS LRTCs play an unequivocal driving role, but rather these dynamics appear to arise through reciprocal interactions between these systems.
Partial correlation analysis shows that neuronal avalanches in distinct T-Δt regimens are independently correlated with either cortical or ANS LRTCs. A, B, The significant values of partial correlations between α(T,Δt) and β(f) with controlling variable βHRV (A) and between α(T,Δt) and βHRV with controlling variable β(f) (B) averaged across multiple frequency bands.
Discussion
Several lines of electrophysiological (Linkenkaer-Hansen et al., 2001), imaging (He, 2011), and behavioral (Palva et al., 2013) evidence show that many features of CNS activity in vivo are scale-free. The scale-free dynamics is relevant because it is a signature characteristic of complex systems poised at criticality (Chialvo, 2010). Operating at a critical state endows the system maximal dynamic range (Kinouchi and Copelli, 2006; Shew et al., 2009), optimal information retention, storage, and transmission capacity (Haldeman and Beggs, 2005; Shew et al., 2011). Nevertheless, several controversies have continued to surround brain criticality and demand further experimental work (Beggs and Timme, 2012; Shew and Plenz, 2013).
The dynamics of a near-critical system can be quantitatively described by power-law scaling exponents (Bak et al., 1987) that in brain dynamics are both predictive of behavioral dynamics (Palva et al., 2013; Smit et al., 2013) and robust biomarkers for many brain diseases (Linkenkaer-Hansen et al., 2005; Nikulin et al., 2012).
However, several distinct forms of scale-free dynamics, including neuronal avalanches in millisecond timescales (Shriki et al., 2013) and LRTCs in 1–100 s scales (Linkenkaer-Hansen et al., 2001), have been observed in human electrophysiological recordings. Both how these phenomena are mutually related and how they interact with similar scale-free dynamics of the ANS (Ivanov et al., 1999) have remained unclear. To address these questions, we set out to obtain a comprehensive quantitative view onto human brain dynamics by using both mesoscopic and macroscopic electrophysiological recordings with SEEG and MEG, respectively. This study advances our understanding of brain criticality in four fronts.
(1) The results indicate that the size distribution of large-amplitude activity cascades (i.e., neuronal avalanches) exhibits robust power-law scaling in both source-reconstructed MEG and white-matter referenced SEEG recordings. These observations are in line with the studies that reported power-law distribution of neuronal avalanche sizes in sensor-level EEG and MEG (Allegrini et al., 2010; Shriki et al., 2013) and invasive ECoG recordings (Solovey et al., 2012; Priesemann et al., 2013). (2) Although there are numerous sensor-level EEG and MEG studies on LRTCs in ongoing human brain activity, no studies have used source-reconstruction methods or SEEG to identify the neuroanatomical characteristics of LRTCs. We show here that the LRTC scaling exponents are systematically distinct both among functional brain systems and between neuronal oscillations in different frequency bands. In most analyses, SEEG and MEG recordings yielded comparable results indicating that neuronal dynamics at spatial meso and macro scales, respectively, are similar. (3) Computational modeling predicts that the exponents of avalanches and LRTCs would be correlated (Poil et al., 2012). We found here a strong correlation between the exponents of neuronal avalanches and LRTCs, which suggests that these two dynamics have shared underlying mechanisms and/or neuronal substrates. (4) We asked whether the critical dynamics of ANS predicts the dynamics in CNS. We found that the scaling exponents of LRTCs in neuronal and HRV time series were correlated in MEG but not in SEEG. This suggests that the impact of ANS dynamics on cortical activity is weak and anatomically widespread so that it becomes observable only in the centimeter-scale coherent activity observable with MEG. Moreover, both CNS and ANS LRTCs had independent partial correlations with neuronal avalanches. These findings together suggest that there is bidirectional interaction between ANS dynamics and macroscopic CNS dynamics.
Data for and against the critical branching process theory for avalanches in neuronal systems
Currently, there are no widely accepted theories for the phenomenology and mechanisms of criticality in neuronal systems. Predictions from the theory of critical branching processes are in line with observations of neuronal dynamics (Beggs and Plenz, 2003), but recent studies also reveal contradictions. In particular, critical branching predicts a constant scaling exponent of 3/2, whereas empirical results indicate that the exponent depends on interelectrode distance (Beggs, 2008), subsampling (Priesemann et al., 2009), and volume conduction (Shriki et al., 2013). Our results show that the scaling exponents of source-reconstructed MEG and SEEG time series are largely above 3/2 for the identical parameters of threshold (T) and time bin width (Δt) that have been used in previous studies (Beggs and Plenz, 2003; Shriki et al., 2013). Taking into account that source reconstructed MEG and SEEG are less affected by signal mixing, these results are in line with the scaling exponent decreasing with increasing linear mixing (Shriki et al., 2013).
The universalism of branching models in the domain of neuroscience is debated because most of the simulation results are mainly consistent with experimental data acquired at submillimeter to millimeter scales (Beggs and Plenz, 2003; Haldeman and Beggs, 2005; Plenz and Thiagarajan, 2007). For instance, the linear relationship between the exponents of avalanche size and lifetime distributions observed at level of individual neurons (Friedman et al., 2012) is in line with earlier studies (Beggs and Plenz, 2004; Plenz and Thiagarajan, 2007). However, such results might not be easily projected onto large-scale neuronal systems, like brains in vivo, where the lifetimes of neuronal avalanches do not exhibit power-law scaling and the size scaling exponents are distinct among neuroanatomical structures (Petermann et al., 2009; Hahn et al., 2010). Furthermore, a growing body of evidence shows that large-scale brain networks are modular and hierarchically organized (Meunier et al., 2010; Gallos et al., 2012). A recent modeling study on avalanches in a modular system showed that avalanche propagation was not confined to the module where it started but still was able to activate only a small subset of neurons in invaded modules (Russo et al., 2014). This suggests that the dynamics of neuronal avalanches in a homogeneous system may be fundamentally different from the dynamics of a system composed of interacting modules.
Neuronal LRTCs are anatomically and spectrally distinct in both MEG and SEEG
Graph theoretical analyses have shown that both the structural connectome (i.e., the predominant axonal pathways) and the dynamically emergent functional connectome of BOLD signal correlations of the human brain are hierarchical and modular networks. Are these internally densely connected modules characterized by dynamics distinct of those in other modules? LRTCs of amplitude of narrowband fluctuations vary across brain areas in EEG and MEG (Linkenkaer-Hansen et al., 2004; Nikulin and Brismar, 2005; Palva et al., 2013), in ECoG (He et al., 2010), and functional systems in fMRI (He, 2011) recordings. To understand the neuroanatomical specificity of LRTCs, we mapped systematically the scaling exponents of LRTCs of source-reconstructed MEG and SEEG time series into a common functional parcellation (Yeo et al., 2011) for the five major frequency bands (δ-, θ-, α-, β-, and γ-band). The largest exponents were found in visual and dorsal attentional networks (α-band) in MEG and SEEG, which is in line with prior EEG (Nikulin and Brismar, 2005) and partially with fMRI (He, 2011) studies. Importantly, SEEG and MEG recordings revealed similar spatio-spectral patterns of the exponents at meso- and macro-spatial-scales, respectively, suggesting that the dynamics reflects a common neuronal process.
However, it is important to note that there are several kinds of noncritical processes that may be associated with power-law LRTCs (e.g., Friedman and Landsberg, 2013; Aitchison et al., 2014; Schwab et al., 2014). Signatures of critical system are also (1) a simple algebraic relationship between the exponents of multiple power laws, (2) the existence of a universal scaling function, and (3) phase transitions caused by tuning of the control parameter (Beggs and Timme, 2012). A modeling study suggests that tuning of the excitation-inhibition balance from overinhibited to underinhibited changes the LRTCs together with a shift in neuronal avalanches from a subcritical to a supercritical regimen, respectively (Poil et al., 2012).
Scaling exponents are correlated at fast- and slow-timescales
Neuronal avalanches and LRTCs coexist in neuronal activity (Beggs and Plenz, 2004; Benayoun et al., 2010; Hahn et al., 2010), but their relationship has remained unclear. We found here that the scaling exponents of avalanche size distribution (α) and LRTCs (β) were correlated in both MEG and SEEG. Surprisingly, the correlation coefficients between α and β were negative for α in the power-law regimen, as tentatively observed previously (Palva et al., 2013), and positive for α in the truncated-power-law regimen. The presence of both negative and positive correlations might be related to LRTCs in the neuronal time series through a simple mechanism. A process with stronger LRTCs facilitates larger avalanches (smaller α), which explain negative correlations at high threshold. At low threshold, a process with weaker LRTCs facilitates larger avalanches (smaller α) mainly made of random events which explain positive correlations accordingly. In time-shuffled surrogate data with no spatial correlations, the positive correlations were largely absent, which is likely associated with the absence of large avalanches caused by random events. In this sense, at meso- and macro-spatial-scales, the LRTCs might play a crucial role in scale-free dynamics of neuronal avalanches.
Footnotes
This work was supported by the Academy of Finland and by the Helsinki University Research Funds. We thank Dr. Hugo Eyherabide for helpful comments on the manuscript.
The authors declare no competing financial interests.
- Correspondence should be addressed to Dr. Alexander Zhigalov, Neuroscience Center, University of Helsinki, P.O. Box 56, Helsinki 00014, Finland. alexander.zhigalov{at}helsinki.fi
