Neuronal avalanches in the resting MEG of the human brain

Abstract

What constitutes normal cortical dynamics in healthy human subjects is a major question in systems neuroscience. Numerous in vitro and in vivo animal studies have shown that ongoing or resting cortical dynamics are characterized by cascades of activity across many spatial scales, termed neuronal avalanches. In experiment and theory, avalanche dynamics are identified by two measures: (1) a power law in the size distribution of activity cascades with an exponent of -3/2 and (2) a branching parameter of the critical value of 1, reflecting balanced propagation of activity at the border of premature termination and potential blowup. Here we analyzed resting-state brain activity recorded using noninvasive magnetoencephalography (MEG) from 124 healthy human subjects and two different MEG facilities using different sensor technologies. We identified large deflections at single MEG sensors and combined them into spatiotemporal cascades on the sensor array using multiple timescales. Cascade size distributions obeyed power laws. For the timescale at which the branching parameter was close to 1, the power law exponent was -3/2. This relationship was robust to scaling and coarse graining of the sensor array. It was absent in phase-shuffled controls with the same power spectrum or empty scanner data. Our results demonstrate that normal cortical activity in healthy human subjects at rest organizes as neuronal avalanches and is well described by a critical branching process. Theory and experiment have shown that such critical, scale-free dynamics optimize information processing. Therefore, our findings imply that the human brain attains an optimal dynamical regime for information processing.

Figure 1.

Discrete MEG events carry correlations of the continuous MEG signal from human resting state. A, Continuous MEG signal of neuronal resting state activity of the human brain (single sensor; NIH). The most extreme point in each excursion beyond a threshold of ±3 SD (horizontal lines) was treated as a discrete event in the signal. Red (green) dots mark positive (negative) events. B, Signal amplitude distributions. The gray curves in the background are the signal amplitude distributions of all individual NIH subjects (based on all channels and all time points). Note that the signal from each sensor was z-normalized by subtracting its mean and dividing by the SD. The blue curve depicts the grand average over all subjects. The red curve depicts the best fit Gaussian distribution for the grand average for the range between ±6 SD. The grand average and the Gaussian fit start deviating from one another around ±2.7 SD. The light blue broken line curve depicts the signal distribution of an empty scanner recording. For clarity, a logarithmic scale is used for the ordinate. The inset depicts the distributions of a single subject and an empty scanner recording using the raw amplitude in pT. C, ETAs for a single sensor. Red/green indicate positive/negative ETA, respectively. D, Discrete events capture most of the significant correlations underlying the continuous MEG signal. Scatter plot shows cross-correlations between original sensor signals from different cortical sites and the corresponding reconstructed signals using ETAs.

Figure 2.

Identification and visualization of spatiotemporal cascades formed by discrete MEG events. A, Raster of events on all sensors (n = 273) in a 10 s segment of recording (single NIH subject). B–D, An example avalanche with cascade size of 20 events, lasting 20 ms and encompassing 19 different sensors. The original time series with identified events (B) leads to the raster of the cascade (C). For visualization, sensors are ordered according to the order of events. A cascade was defined as a series of time bins in which at least one event occurred, ending with a silent time bin. Here the time bin width was 3.3 ms, twice the sampling time step (1.67 ms; 600 Hz). This cascade is visible as a positive-signed propagation in the lower left part of the sensor array (D). In each panel of D, the black dots mark which sensors were active in that time bin. The last panel depicts the set of all sensors that participated in the cascade.

Figure 3.

Cascade size distributions follow power laws, as expected for neuronal avalanches. A–C Cascade size distributions for a single NIH subject using axial gradiometers and Δt = 3.3 ms. A, Solid black line depicts original data and dashed red line corresponds to phase-shuffled data. Dashed black line represents a power law with an exponent of −3/2. Arrow indicates the number of sensors (N) in the analysis (system size). B, Cascade size distributions for subsamples of the sensor array. Line color and arrows indicate the number of sensors (N) in the analysis. Upper right inset: Diagrams of the sensor array with colored subsamples. Lower left inset: The same cascade size distributions, plotted as a function of the scaled axis Z = S/N. C, Cascade size distributions for the coarse-grained array. Black line indicates original data; green line, coarse-grained data. Inset: Diagram of coarse-grained sensor array with sensors grouped in clusters of ∼4 sensors. D–F, as in A–C for the same NIH subject but using virtual planar sensors and Δt = 3.3 ms. G–I, as in A–C for a Cambridge subject using planar sensors and Δt = 12 ms.

Figure 4.

Cascade size distributions of empty scanner data do not display power law behavior. Cascade-size distributions from 8 empty scanner recordings (at NIH) segmented with Δt = 3.3 ms and a threshold of ±3 SD. The distributions fall off much earlier than the size of the array (n = 273 sensors).

Figure 5.

Neuronal avalanches in human MEG reveal power law exponent of α = −3/2 at critical branching parameter σ = 1. A, Phase plots of the power law exponent, α, versus the branching parameter, σ, using axial sensors. Each point represents a single subject at a single Δt, where different colors correspond to different values of Δt (see color key; n = 104 subjects). B–D, Average phase plots of the exponent, α, versus the branching parameter, σ, for the several array types examined. Vertical and horizontal bars denote SD. Solid vertical and horizontal lines denote the point σ = 1, α = −3/2. Insets depict the corresponding sensor arrays. B, Axial (black circles) and virtual planar (blue squares) sensors. C, Subsamples of the array (only for axial sensors). Error bars were omitted for clarity of presentation. D, Coarse-grained array for axial (black circles) and planar (blue squares) sensors. E, F, Robustness to changes in threshold and peak detection method. E, Average phase plots from all sensors for threshold values ranging from ±2.7 to ±4.2 SD. Here and throughout the manuscript, an event was identified as the most extreme point in each excursion beyond the threshold (inset). Increasing the threshold from ±2.7 to ±4.2 SD reduced the avalanche rate by an order of magnitude from 27.3 to 2.4 Hz for bin width of 3.3 ms. F, Same as in E, but with a peak detection method that identifies events at all local extremum points beyond the threshold (inset). The change in peak detection method did not change substantially the overall rate of avalanches, which was 27.1 Hz at a threshold of ±2.7 SD and 2.5 Hz at ±4.2 SD.

Figure 6.

Behavior of power law exponent and branching parameter as a function of timescale Δt for Cambridge subjects. A–D, As in Figure 5 but based on MEG recordings at the Cambridge facility. A, Phase plots of α vs σ. Each point represents a single subject at a single Δt. Different colors correspond to different values of Δt (see color key). B–D, Average phase plots of α versus σ for the several array types examined. Vertical and horizontal bars denote SD. Solid vertical and horizontal lines denote the point σ = 1, α = −3/2. Insets depict the corresponding sensor arrays. B, Original (planar) sensors. C, Subsamples of the array. Error bars were omitted for clarity of presentation. D, Coarse-grained array. E, Avalanche statistics of Cambridge subjects with eyes closed, original recordings (black, circles) and with eyes closed, with the independent component most associated with eye movement removed (blue, squares). F, Avalanche statistics of Cambridge subjects with eyes closed (black, circles) and open (blue, squares). In both cases, the independent components most associated with eye movement were removed.

Figure 7.

Robustness of power law exponent and branching parameter over time and between recording sessions. A, B, Correlation of the power law exponent α (A) and branching parameter σ (B) for NIH subjects (n = 104) comparing the first 2 min of each 4 min recording with the last 2 min. C, D, as in A, B for Cambridge subjects (n = 20) comparing two visits separated by at least 1 week.

Figure 8.

Neural network simulations to study the effect of linear mixing at the sensors on estimating avalanche dynamics. A, Illustration of the neuron grid and the sensor array. The distance between neighboring sensors, d, was four times the distance between neighboring neurons. The full network had N = 2025 neurons arranged in a 45 × 45 grid, M = 64 sensors arranged in an 8 × 8 grid, and the margins on each side of the sensor array were 8 times the distance between neighboring neurons (see Materials and Methods). B, One-dimensional illustration of the overlap between two adjacent sensors on the two-dimensional array. Each sensor had a Gaussian weight profile and the SD was varied to explore the effect of small (SD = 0.5 d) versus large (SD = 1.5 d) overlap. C, Cascade size distribution from a simulation of independent Poisson processes on the neuron grid. Dashed black line represents a power law with an exponent of −3/2. D, Cascade size distributions from the simulation in C as captured by the sensor array for different levels of sensor overlap. E, Cascade size distribution from a simulation of the network at criticality as observed on the neuron grid. F, Cascade size distributions from the simulation in E as captured by the sensor array for different levels of sensor overlap. G–H, Dependence of the power law exponent (G) and the branching parameter (H) on the level of sensor overlap.

Figure 9.

PLI analysis yields similar results for human and empty scanner data. A, B, PLI distributions of a single human subject (A) and a single empty scanner (B), both at the Cambridge MEG facility. The recordings were filtered using Hilbert wavelet pairs at four wavelet scales, with corresponding frequency bands included in the legend. C–H, PLI analysis using band-pass filtering and the Hilbert transform to acquire phase information. C, D, Single human and empty scanner analysis (Cambridge). E, F, Single human and empty scanner analysis (NIH). G, H, Averages of 104 human NIH recordings (G) and 8 NIH empty scanner recordings (H). SEM was narrower than line width.