Statistical modeling of adaptive neural networks explains co-existence of avalanches and oscillations in resting human brain

Neurons in the brain are wired into adaptive networks that exhibit collective dynamics as diverse as scale-specific oscillations and scale-free neuronal avalanches. Although existing models account for oscillations and avalanches separately, they typically do not explain both phenomena, are too complex to analyze analytically or intractable to infer from data rigorously. Here we propose a feedback-driven Ising-like class of neural networks that captures avalanches and oscillations simultaneously and quantitatively. In the simplest yet fully microscopic model version, we can analytically compute the phase diagram and make direct contact with human brain resting-state activity recordings via tractable inference of the model’s two essential parameters. The inferred model quantitatively captures the dynamics over a broad range of scales, from single sensor oscillations to collective behaviors of extreme events and neuronal avalanches. Importantly, the inferred parameters indicate that the co-existence of scale-specific (oscillations) and scale-free (avalanches) dynamics occurs close to a non-equilibrium critical point at the onset of self-sustained oscillations.


I. INTRODUCTION
Synchronization is a key organizing principle that leads to the emergence of coherent macroscopic behaviors across diverse biological networks [1].From Hebb's "neural assemblies" [2] to synfire chains [3,4], synchronization has also strongly shaped our understanding of brain dynamics and function [5].The classic and arguably most prominent example of large scale neural synchronization are brain oscillations, first reported about a century ago [6]: periodic, large deflections in electrophysiological recordings, such as electroencephalography (EEG), magnetoencephalography (MEG), or local field potential (LFP) [6,7].Because oscillations are thought to play a fundamental role in brain function, their mechanistic origins have been the subject of intense research.According to the current view, the canonical circuit that generates prominent brain rhythms such as the alpha oscillations and the alternation of up-and down-states utilizes mutual coupling between excitatory (E) and inhibitory (I) neurons [8][9][10].Alternative circuits, including I-I population coupling, have been proposed to explain other brain rhythms such as high-frequency gamma oscillations [11][12][13].Setting biological details aside, the majority of research has predominantly focused on the emergence of synchronization at a preferred temporal scale -the oscillation frequency.
Brain activity also exhibits complex, large-scale cooperative dynamics with characteristics that are antithetic to those of oscillations.In particular, empirical observations of "neuronal avalanches" have shown that brain rhythms coexist with activity cascades in which neuronal groups fire in patterns with no characteristic time or spa-tial scale, suggesting that the brain may operate near criticality [14][15][16][17][18][19][20][21][22][23][24][25].In this context, the coexistence of scale-free neuronal avalanches with scale-specific oscillations suggests an intriguing dichotomy that is currently not understood.On the one hand, models of brain oscillations are very specific and seek to capture physiological mechanisms underlying particular brain rhythms.On the other hand, attempts to explain the emergence of neuronal avalanches almost exclusively focus on criticalityrelated aspects and ignore the coexisting behaviors such as oscillations, even though they themselves may be constitutive for understanding the putative criticality.Among the few exceptions [26][27][28][29][30], Poil et al proposed a probabilistic integrate and fire (IF) spiking model with E and I neurons that generates long-range correlated fluctuations reminiscent of MEG oscillations in the resting state, with supra-threshold activity following power-law statistics consistent with neuronal avalanches and criticality [26].More recently, by adopting a coarse-grained Landau-Ginzburg approach to neural network dynamics, Di Santo et al have shown that neuronal avalanches and related putative signatures of criticality co-occur at a synchronization phase transition, where collective oscillations may also emerge [28].These results were successively extended to a hybrid-type synchronization transition in a generalized Kuramoto model [31].
While both these and other proposed approaches show that neuronal avalanches may coexist with some form of network oscillations [26,30] or network synchronization [28,31], they suffer from three major shortcomings.First, these models are neither simple (e.g., in terms of parameters) nor analytically tractable, making an exhaustive exploration of their phase diagram out of reach.
Second, none of the two models simultaneously captures events at the microscopic scale (individual spikes) and macroscopic scale (collective variables).Third, it is not clear how to connect these models to data rigorously, beyond relying on qualitative correspondences.
Here we propose a minimal, microscopic, and analytically tractable model class that can capture a wide spectrum of emergent phenomena in brain dynamics, including neural oscillations, extreme event statistics, and scale-free neuronal avalanches [14].This model class is inspired by the recent theoretical observation that systems with many interacting degrees of freedom and coexistence of distinct phases may develop self-oscillations in the presence of feedback loops between control and order parameters [32].Brain dynamics, with its patterns of state transitions [33,34], long-range correlated fluctuations and cooperative behaviors showing marks of criticality [16,[35][36][37][38], relies on a number of regulatory feedback loops that are considered at the core of its function [10,39].Here, we hypothesize that basic feedback mechanisms controlling the excitability of the system could produce the coexistence of the antithetic scalespecific oscillations and scale-free avalanches in brain activity.To test this hypothesis, we build upon the wellestablished analogy between neural networks and spin systems [38,[40][41][42][43][44], and put forward a non-equilibrium extensions of the Ising model of statistical physics with an extra feedback loop which enables self-adaptation.As a consequence of feedback, neuronal dynamics is driven by the ongoing network activity, generating a rich repertoire of dynamical behaviors.As in previous applications of the static Ising model to neural systems [42][43][44], spins model the spiking, binary nature of individual neurons (firing versus silent).Thus, the structure of the simplest model from this class permits microscopic network dynamics investigations as well as analytical mean-field solution in the Laudau-Ginzburg spirit, and in particular, allows us to construct the model's phase diagram.
The tractability of our model enables us to make direct contact with MEG data on the resting state activity of the human brain.With its two free parameters inferred from data, the model closely captures brain dynamics across scales, from single sensor MEG signals to collective behavior of extreme events and neuronal avalanches.Remarkably, the inferred parameters indicate that scale-specific (neural oscillations) and scale-free (neuronal avalanches) dynamics in brain activity coexist close to a non-equilibrium critical point that we proceed to characterize in detail.

A. Adaptive Ising model
We consider a population of interacting neurons whose dynamics is self-regulated by a time-varying field that depends on the ongoing population activity level (Fig. 1A).The N spins s i = ±1 (i = 1, 2, ..., N , N = 10 4 in our simulations unless specified differently) represent excitatory neurons that are active when s i = +1 or inactive when In the simplest, fully homogeneous scenario described here, neurons interact with each other through synapses of equal strength J ij = J = 1.The ongoing network activity is defined as , as the magnetization of the Ising model) and each neuron experiences a uniform negative feedback h that depends on the network activity as ḣ = −cm, with c determining the strength of the feedback.Neurons s i are stochastically activated according to the Glauber dynamics, where the new state of neuron s i is drawn from the marginal Boltzmann-Gibbs distribution P (s i ) ∝ exp(β hi s i ), with hi = j =i J ij s j +h, where β is reminiscent of the inverse temperature for an Ising model (see Appendix B).
Multiple interpretations of this model are possible.On the one hand, negative feedback can be identified with a mean-field approximation to the inhibitory neuron population that uniformly affects all excitatory neurons with a delay given by the characteristic time c −1 (see Appendix B).On the other hand, feedback could be seen as intrinsic to excitatory neurons, mimicking, e.g., spike-threshold adaptation [45][46][47][48].Exploration-worthy (and possibly more realistic) extensions within the same model class are accessible by considering two ways in which geometry can enter the model.First, as in the standard Ising magnet, the interactions J can be restricted to simulate local excitatory connectivity, e.g., to nearest neighbors on a 2D lattice.Second, feedback h i to neuron i could be derived from a local magnetization in a neighborhood around neuron i instead of the global magnetization; in the interesting limiting case where ḣi = −cs i , each neuron would feed back on its own past spiking history only, and the model would reduce to a set of coupled "binary oscillators" [32].Irrespective of the exact setting, the model's mathematical attractiveness stems from its tractable interpolation between stochastic (spiking of excitatory units) and deterministic (feedback) elements.
Network behavior is determined by feedback strength c and inverse temperature β.In the fully-connected continuous-time limit, the model can be described with the following Langevin equations: where ξ is unit uncorrelated Gaussian noise; the stochastic term thus has amplitude b = 2/(βN ).Equations (1) can be linearized around the stationary point (m * = 0, h * = 0) to calculate dynamical eigenvalues and construct a phase diagram (Fig. 1B): For c = 0, h = 0, the model reduces to the standard infinite-dimensional (mean field) Ising model with a second order phase transition at β = β c = 1.At non-zero feedback, c > 0, the model is driven out of equilibrium and its critical point at β c coincides with an Andronov-Hopf bifurcation [32,49].For c below a threshold value c * = (β − 1) 2 /4β, m(t) is described by an Ornstein-Uhlenbeck process (O-U) independently of β.For β < β c , the system is stable and shows a crossover from a stable node with exponential relaxation (two negative real eigenvalues) to a stable focus with oscillation-modulated exponential relaxation (two imaginary eigenvalues; "resonant regime") when c increases beyond c * (Fig. S1).In the resonant regime, c > c * , oscillations become more prominent as the critical point β c = 1 is approached, finally transitioning into self-sustained oscillations for β > β c (Fig. S2).
We focus on the resonant regime below and at the critical point, and study the reversal times and zero-crossing areas of the total network activity m(t) (Fig. 1C).The distribution P (a 0 ) of the zero-crossing area follows a power-law behavior with an exponent τ = 1.29 ± 0.01 in the vicinity of the critical point.As β decreases, the scaling regime shrinks until it eventually vanishes for small enough β.Similar behavior is observed for the distribution P (t) of reversal times.This distribution also follows a power-law with an exponent α t = 1.40 ± 0.01 near the critical point (Fig. 1D).Both distributions have an exponential cutoff related to the characteristic time of the network activity oscillations, 1/c; this cutoff transforms into a hump as β → 1 and c c * (β), i.e., as oscillations in m(t) become increasingly prominent (Fig. S3).Importantly, for the non-interacting (J = 0) model, the distributions P (a 0 ) and P (t) follow a purely exponential behavior (Fig. 1D, inset), indicating that the coexistence of oscillatory bursts and power-law distributions for the network activity requires neuron interactions as well as the adaptive feedback (Fig. S4).

B. Model inference from local resting-state brain dynamics
In the resonant regime below the critical point (c > c * , β < β c ), it is possible to analytically compute the autocorrelation function of the ongoing network activity m(t) in the linear approximation [50]: where γ = (1 − β)/2 and ω = βc − (1 − β) 2 /4.The autocorrelation C(τ ) can be used to infer model parameters β and c from empirical data by moment matching, thereby locating the observed system in the phase diagram (Fig. 1B).
We test the proposed approach on MEG recordings of the awake resting-state of the human brain (Appendix A).We first analyze brain activity on individual MEG sensors.To this end, we compare the magnetic field recorded on individual MEG sensors with the magnetization m of the model (Fig. 1).This analogy relies on the nature of the brain magnetic fields captured by the MEG, which are generated by synchronous postsynaptic currents in cortical neurons [51].Because the intensity of such currents depends on the number of actively firing neurons, the temporal fluctuations in the MEG signal are related to average firing rate.Similarly, in the model, the fluctuations in the magnetization m are directly related to the number of active neurons over time, that is the firing rate.
During resting wakefulness the brain activity is largely dominated by oscillations in the alpha band (8 − 13 Hz) (Fig. 2A), which has been the starting point of many investigations [7,34,52,53] including ours reported below; similar results are also obtained for the broadband activity (Fig. S5).After isolating the alpha band, we estimate the quantities γ and ω by fitting the empirical C(τ ) to the functional form given by Eq (3).Fig. 2B illustrates the typical quality of the fit and the qualitative resemblance between the model and MEG sensor signal dynamics.
Since our model is fit to reproduce the second-order statistical structure in the signal, we next turn our attention to signal excursions over threshold, a higher-order statistical feature routinely used to characterize bursting brain dynamics [18,34,37,[54][55][56].To that end, we construct the distribution of (log) areas under the signal above a threshold ±e (Fig. 2C) [34].P (log a e ) is bell-shaped, featuring strongly asymmetric tails for MEG sensors as well as the model (Fig. 2C).Variability across subjects is mostly related to signal amplitude modulation, resulting in small horizontal shifts in P (log a e ) but no variability in the distribution shape.Remarkably, the rescaled distribution is independent of the threshold e over a robust range of values, and is well-described by a Weibull form, 2C, bottom panel inset; Fig. S6).Taken together, these observations indicate that our model has the ability to capture non-trivial aspects of amplitude statistics in MEG signals, within and across different subjects (Fig. S7).
Parameters inferred across all sensors and subjects sug-gest baseline values of β = 0.99 and c = 0.01 that are well matched to data, which we use for all subsequent analyses (unless stated otherwise).Specifically, we find the best-fit β values strongly concentrate in a narrow range around β ≈ 0.99 (β = 0.986 ± 0.006; c = 0.012 ± 0.001), very close to the critical point (Fig. 2D and Fig. S8).Even though all analyzed signals are bandpass-limited to a central frequency around 10 Hz by filtering, closeness to criticality appears to strongly correlate with the fraction of total power in the raw signal in the alpha band (R 2 = 0.59; p < 0.001).This suggests that alpha oscillations may be closely related to critical brain tuning during the resting state [19,35,36,52,55].
A classic fingerprint of tuning to criticality is the emergence of long-range temporal correlations (LRTC), which have been documented empirically [35-37, 52, 55, 57, 58].LRTC in the alpha band have been investigated primarily by applying the detrended fluctuations analysis (DFA) to the amplitude envelope of MEG or EEG signals in the alpha band (Appendix C) [26,35,55,59].Briefly, DFA estimates the scaling exponent α of the root mean square fluctuation function F in non-stationary signals with polynomial trends [60].For signals exhibiting positive (or negative) LRTC, F scales as F ∝ n α with 0.5 < α < 1 (or 0 < α < 0.5, respectively); α = 0.5 indicates the absence of long range correlations; α also approaches unity for a number of known model systems as they are tuned to criticality [61,62].
To test for the presence of LRTC using DFA, we analyzed the scaling behavior of fluctuations and extracted their scaling exponent α [35,52].To avoid spurious correlations introduced by signal filtering, α was estimated over the range 2 s < n < 60 s (Fig. 2E) [35,52].We find that α is consistently between 0.5 and 1 for all MEG sensors and subjects, in agreement with previous analyses [35,55,57,58,63,64].Importantly, model-free α values measured across MEG sensors positively correlate with the inferred β values from the model (Fig. 2F), indicating that higher β values are diagnostic about the presence of long-range temporal correlations in the amplitude envelope.
Taken together, our analyses so far show that the adaptive Ising model recapitulates single MEG sensor dynamics by matching their autocorrelation function and the distribution of amplitude fluctuations, and further suggest that the true MEG signals are best reproduced when the adaptive Ising model is tuned close to, but slightly below, its critical point (β 1).

C. Scale invariant collective dynamics of extreme events
We now turn our attention to phenomena that are intrinsically collective: (i) coordinated supra-threshold bursts of activity, which emerge jointly with LRTC in alpha oscillations [26,55]; and (ii) neuronal avalanches, i.e., spatio-temporal cascades of threshold-crossing sen-sor activity, which have been identified in the MEG of the resting state of the human brain [19,36].Both of these phenomena are generally seen as chains of extreme events that are diagnostic about the underlying brain dynamics [18,54,[65][66][67][68].
We start by defining the instantaneous network excitation, A (t), as the number of extreme events co-occurring within time bins of size across the entire MEG sensor array (Appendix D).For each sensor, extreme events are the extreme points in that sensor's signal that exceed a set threshold e (Fig. 3A).For a given threshold, network excitation A depends on the size of the time bin that we use to analyze the data (Fig. 3B).To make contact with the model, we parcel our simulated network into K equally-sized disjoint subsystems of n sub = N/K neurons each, and consider each subsystem activity m µ , µ = 1, . . ., K, as the equivalent of a single MEG sensor signal.Network excitation, A , for the model then follows the same definition as for the data, allowing us to perform direct side-by-side comparisons of extreme event statistics.
We first study the distribution of network excitation, P (A ).We set e = 2.9SD both for MEG data and for the model [19].Even though P (A ) generally depends on , the distributions corresponding to different collapse onto a single, non-exponential master curve when A is rescaled by A , the average instantaneous network excitation (Fig. 3C).Excitation distribution is thus invariant under temporal coarse-graining and the number of extreme events scales non-trivially with , in contrast to phase-shuffled surrogate data (Appendix E).Remarkably, model simulations fully recapitulate this data collapse as well as the non-exponential extreme event statistics.
Model simulations reproduce the distribution of network excitation, P (A ), to within the variability observed among subjects (Fig. 3C, inset), for given values of .To quantify how close the model distribution P m (A ) is to the data-derived average distribution P d (A ), we calculate the Kullback-Leibler (KL) divergence [69]: This is to be compared with the average KL divergence across subjects: , averaged across all pairs of MEG subjects indexed by i and j.The data-model divergence is very small and within the range of variability across subjects (D dm (A ) D dd (A ), see Table 1), suggesting that the model quantitatively reproduces the measured distributions to the degree that can be expected given natural variability in the data.
Periods of excitation (A = 0) are separated by periods of quiescence (A = 0) of duration I = n , where n is the number of consecutive time bins with A = 0.The distribution of quiescence durations, P (I ), is in- MEG Model Surrog  variant under temporal coarse-graining when rescaled by I , the average quiescence duration, collapsing onto a single, non-exponential master curve (Fig. 3D).As was the case with the distribution of network excitation, the model-predicted distribution of quiescence durations, P m (I ), also diverges from the data average P d (I ) by an amount that is within the range of variability among subjects (Fig. 3D, inset; Table 1).
We furthermore show that the overall probability P 0 ( ) of finding a quiescent time bin follows a non-exponential relation P 0 ( ) = exp −a β I , with β I 0.6 (Fig. 3D, lower inset), indicating that extreme events grouped into bins of increasing size are not independent [70].These results are robust to changes in N and n sub so long as the number of subsystems K is fixed or does not change considerably (Figs S9-S10); otherwise, the threshold e that defines an extreme event should be adjusted accordingly (Fig. S11), in particular, to closely reproduce the distribution of quiescence durations P (I ) (Fig. S11).Finally we notice that the quantities A and I scale as a power of the bin size (Fig. S12), and are connected to each other by a relationship of the form A ∼ I b AI (Fig. S12).This implies that, for a fixed value of e, both the distribution P (A ) and P (I ) are controlled by a sin-gle quantity, e.g., the average network excitation A .
In sum, our simple model at baseline parameters provides a robust account of the collective statistics of extreme events.We emphasize that the excellent match to the observed long-tailed distributions is only observed for the inferred value β 0.99 very close to criticality; already for β = 0.98, we observe significant deviations from data (Figs.S13, S14), demonstrating that excitation and quiescence distributions represent a powerful benchmark for collective brain activity.

D. Scale-free neuronal avalanches occur concomitantly with oscillations
A neuronal avalanche is a maximal contiguous sequence of time bins populated with at least one extreme event per bin (Fig. 4A) [14,19]; every avalanche thus starts after and ends with a quiescent time bin (A = 0) (see Appendix D for details).Typically, neuronal avalanches are characterized by their size s, defined as the total number of extreme events within the avalanche.Avalanche sizes have been reported to have a scale-free power-law distribution [14,16,19,24,36].
We estimate the distribution P (s) of avalanche sizes in the resting state MEG, and compare it with the distribution obtained from model simulation at close-tocritical baseline parameter set (Fig. 4B).Both distributions are described by a power-law with an exponential cutoff [19] and show an excellent match across subjects and for individual subjects.Again, the KL divergence between the mean empirical and the model distribution is smaller than the mean KL divergence estimated among MEG subjects (Table 1).Phase scrambled surrogate data strongly deviate from the power-law observations, as do model predictions when parameter β is moved even marginally below 0.99 (Fig. S15).These results are independent of the N and n sub so long as the number of subsystems K is fixed (Fig. S16) or does not change considerably (Fig. S11).Importantly, the model also reproduces the scaling relation s (d) ∼ d ζ that connects average avalanche sizes s and durations d.Unlike the power-law exponent of avalanche size distribution that typically depends on time bin size [14,36], the exponent ζ does not depend on , as shown by the data collapse for both MEG data and model (Fig. 4B, the inset).While the scaling behavior is reproduced qualitatively, the inferred and model-derived values of ζ are not in quantitative agreement, likely due to the overly simplified mean-field connectivity assumed by our model.

III. DISCUSSION
In this paper, we put forward the adaptive Ising class of models for capturing large scale brain dynamics.Quite generally, these models combine a microscopic and stochastic description of excitatory neuron spiking with coarse-grained mean-field model of activity-dependent feedback.This endows the models with several unique characteristics that we discuss in turn: (i) the ability to generate a diverse range of stylized behaviors observed in real brain dynamics and locate these behaviors in the model's phase diagram; (ii) the possibility to connect the model to a wide range of known theoretical results in statistical physics and theoretical neuroscience; (iii) the ability to rigorously infer model parameters from data, quantitatively test its predictions across a range of spatial and temporal scales, and thus derive biological insight about brain function.4)).

A. Diversity of dynamical behaviors in a simple non-equilibrium model
By combining local interactions with a time-varying field in the form of an activity-dependent feedback, the adaptive Ising model exhibits a phase transition to a selfoscillatory behavior at the Andronov-Hopf bifurcation critical point that inherits the characteristics of the classical Ising ferromagnet second-order phase transition [32].In proximity and slightly below the critical point, the ongoing network activity m -the equivalent of the Ising magnetization -shows intermittent oscillations, whose associated reversal time and zero crossing area are powerlaw distributed.To our knowledge, this is the simplest model class that reproduces the stylized coexistence of neuronal avalanches and oscillations, the two antithetic features of real brain dynamics.Moreover, these features coexist already in the most basic, mean-field formulation of the model, whose phase diagram can be computed analytically.Interestingly, in this formulation, individual units are neither intrinsic oscillators themselves [31,39], nor are they mesoscopic units operating close to a Hopf bifurcation [10,39,71], and the collective dynamics is therefore not a result of oscillator synchronization (even though this regime can be captured as well by a different realisation of an adaptive Ising model).Our proposal thus provides an analytically-tractable alternative toor perhaps a reformulation of-existing models [26,28], which typically implicate particular excitation/inhibition or network resource balance to open up the regime where oscillations and avalanches may coexist.

B. Connections to results in statistical physics
Starting with the seminal work of Hopfield [41], the functional aspects of neural networks have traditionally been studied with microscopic spin models or attractor neural networks.These systems qualitatively reproduce the associative memory functionality ascribed to real neural networks [72] and their phase transitions have been thoroughly studied with statistical mechanics tools [73].The associated inverse (maximum entropy) problem recently attracted great attention in connecting spin models to data [42,74], in particular with regard to criticality signatures [38].While the initial formulations of the maximum-entropy problem aimed at modeling the stationary distribution from which activity patterns were drawn independently, later work specifically focused on the temporal correlation structure in the neural activity, by modeling across-time interactions between individual neurons [75][76][77].Parallel work meanwhile considered spin-models responding to global known modulation, e.g., via stimulus, in the stimulus-dependent maximum-entropy approach [78], or global latent modulation which can result in critical (Zipf-like) stationary distributions of activity [79][80][81], or dynamics reminiscent of avalanches [82].
Despite such generalizations, the dynamical expressive power of maximum-entropy stationary, kinetic, or latentvariable models has been limited: while detailed patterns of short-term temporal correlations and statistical criticality could be captured from data, the rhythmic behavior of brain oscillations was beyond the practical scope of these models.Adaptive Ising model class can be seen as a natural, yet orthogonal, extension to the previous work that enables oscillations and furthermore permits us to explore an interesting interplay of mechanisms, for example, by having self-feedback drive Hopfield-like networks (with memories encoded in the coupling matrix J) through sequences of stable states.
Another link to rich existing literature becomes apparent when we consider connectivity degrees of freedom (encoded in J and in the pooling that drives the feedback; Appendix B) in order to model cell types (e.g., inhibitory vs excitatory), the spatial structure of the cortex, or to capture empirically established topological features of real neural networks.Within equilibrium statistical mechanics, the effects of various lattices or disordered connectivity in spin systems have been studied rigorously; phenomenological arguments, supported by preliminary results, suggest that adding the feedback that drives the adaptive Ising model out of equilibrium does not change the fundamental nature (e.g., critical exponents) of the critical point.This permits us to directly translate existing equilibrium results into the non-equilibrium setup: an example would be the introduction of scale-free topology among excitatory neurons into our model (example in the Appendix B).Looking towards the future, powerful tools of statistical physics, such as the Landau-Ginzburg theory and the renormalization group, can be brought to bear on the generalizations of the non-equilibrium adaptive Ising model, enabling rapid progress based on existing equilibrium results.

C. Inferring the model to derive biological insights
Since the adaptive Ising model can be solved analytically in the mean field limit, we can infer its two key parameters by matching the auto-correlation function of MEG signals to the model-predicted auto-correlation function.In contrast to previous work [26,28], we do not make contact with existing data by qualitatively matching the phenomenology, but by proper parameter inference.The inferred parameters consistently place the model very close to its critical point, supporting the hypothesis that alpha oscillations represent brain tuning to criticality [10,34,35,52,55].The possibility of mapping empirical data to a defined region in the adaptive Ising model phase diagram through parameter inference paves the way for further quantification of the relationship between measures of brain criticality and healthy, developing, or pathological brain dynamics along the lines developed recently [64,66,83].
Our inferred model makes a wide range of further predictions that can be confronted with data.By parcelling the simulated system into groups of spins, we can mimic the signals captured by multiple MEG sensors over the cortex.We find a remarkable quantitative correspondence between the non-exponential and scale-invariant distribution of network excitation observed across the MEG sensors and in our simulation, which holds when averaged across individuals or within single individuals.Our analysis clearly demonstrates that extreme events are non-independent across space and time.This nontrivial spatio-temporal organization of extreme events is strongly indicative of a network state close to criticality [18,65]; the agreement between model and data breaks down for surrogate data as well as for models further removed from criticality.Moreover, the extreme events coalesce into neuronal avalanches as reported previously [14,16,19,36], which our model reproduces as well.Taken together, our model provides a remarkably broad account of brain dynamics across spatial and temporal scales.
Despite these successes, we openly acknowledge the quantitative failures of our model: (i) at the single sensor level, small deviations exist in the distributions of log activity (Fig. 2C), likely due to very long timescales or non-stationarities in the MEG signals [19,84,85]; (ii) small deviations beyond the range of data variability exist in the probability distribution P (I ) for a narrow range of intermediate quiescence intervals, even when the rest of the distribution is reproduced very well (Fig. 3D, inset); (iii) the scaling exponent governing the relation between the avalanche size and duration, ζ, is not reproduced quantitatively (Fig. 4B, inset).Furthermore, beyond the two key model parameters that were inferred directly from individual sensors (β, c), quantitative data analysis of extreme events requires additional parametric choices (time bin , threshold e, system size N and subsystem size n sub ), both for empirical data as well as model simulations.While we successfully demonstrate the scaling invariance of the relevant distributions with respect to and robustness with respect to N and n sub at fixed N/n sub , a close match to data still requires choosing one extra parameter (e.g., threshold e).Concerning this point, we have verified that an even closer agreement between empirical and numerical distributions can be achieved setting slightly different threshold values e on MEG and model data-in particular, when the thresholds are adjusted so that the data and the model are perfectly matched in the average activity rate.
Despite these valid points of concern, we find it remarkable that such a simple and tractable model can quantitatively account for so much of the observed phenomenology.Future work should first consider connectivity beyond the simple all-to-all mean-field version that we introduced here, likely leading to a better data fit and new types of dynamics, e.g, cortical waves.Preliminary model simulations show that the exponent ζ, which connects avalanche sizes and durations, is affected by the connectivity and, furthermore, more closely matches the value of ∼ 1.3 characteristic of the data for nearestneighbor 2D connectivity on a square lattice.Second, we strongly advocate for rigorous and transparent data analysis and quantitative-not only stylized-comparisons to data.To this end, care must be taken not only when inferring the essential model parameters beyond the linear approximation [86], but also when treating the hidden "degrees of freedom" related to data analysis (specifically, subsampling, temporal discretization, thresholding etc.) [14,19,36,87,88].Third, it is important to confront the model with different types of brain recordings; a real success in this vein would be to account simultaneously for the activity statistics at the microscale (spiking of individual neurons) as well as at the mesoscale (coarsegrained activity probed with MEG, EEG, or LFP).
with NIH guidelines for human subjects.The sampling rate was 600 Hz, and the data were band-pass filtered between 1 and 150 Hz.Power-line interferences were removed using a 60 Hz notch filter designed in Matlab (Mathworks).The sensor array consisted of 275 axial first-order gradiometers.Two dysfunctional sensors were removed, leaving 273 sensors in the analysis.Analysis was performed directly on the axial gradiometer waveforms.The data analyzed here were selected from a set of MEG recordings for a previously published study [19], where further details can be found.For the present analyses we used the subjects showing the highest percentage of spectral power in the alpha band (8)(9)(10)(11)(12)(13).Similar results were obtained for randomly selected subjects.

APPENDIX B: FURTHER DETAILS ON THE ADAPTIVE ISING MODEL
The model is composed of a collection of N spins s i = ±1 (i = 1, 2, ..., N ) that interact with each other with a coupling strength J ij .In our analysis, the N spins represent excitatory neurons that are active when s i = +1 or inactive when s i = −1, and J ij > 0. Furthermore, we consider the fully homogeneous scenario, with neurons interacting with each other through synapses of equal strength J ij = J = 1.However, interesting generalization with non-homogeneous, negative, nonsymmetric J ij are possible, to include, for example, the effect of inhibitory neuronal population and structural and functional heterogeneity.The s i are stochastically activated according to the Glauber dynamics, where the state of a neuron is drawn from the marginal Boltzmann-Gibbs distribution The spins experience an external field h, a negative feedback that depends on network activity according to the following equation, where c is a constant that controls the feedback strength, and the sum runs over a neighborhood of the neuron i specified by N i ; index j enumerates over all the elements of this neighborhood.Depending on the choice of N i , the feedback may depend on the activity of the neuron i itself (self-feedback), its nearest neighbors, or the entire network -the case we considered in the main paper.In a more realistic setting including both excitatory (J ij > 0) and inhibitory neurons (J ij < 0), one could then take into account the different structural and functional properties of excitatory and inhibitory neurons by considering different interaction and feedback properties [89].In our simulations, one time step corresponds to one system sweep -i.e.N spin flips -of Monte Carlo updates, and Eq (6) is integrated using ∆t = 1/N .Note that this choice of timescales for deterministic vs stochastic dynamic is important, as it interpolates between the quasiequilibrium regime where spins fully equilibrate with respect to the field h, and the regime where the field is updated by feedback after each spin-flip and so spins can constantly remain out of equilibrium.∆t is generally much smaller than the characteristic time of the adaptive feedback that is controlled by the parameter c.

Mapping between the adaptive Ising model and an E-I network
We first encode a classic E-I model that leads to sustained oscillations, into an Ising model framework.We can think of this model as a single network of N units whose coupling matrix J ij is asymmetric and is structured into two blocks that correspond to an excitatory and an inhibitory subpopulation.Specificcally, the network consists of a population 1, which is self-exciting with strength J 11 and which excites population 2 with strength J 12 , while population 2 is inhibiting population 1 with strength J 21 .
It can be demonstrated that the mean-field dynamics of this stochastic system of Ising-like neurons in the limit of large populations is described by a Liouville deterministic equation of the form (m i i = 1, 2 is the average spiking rate of population i) [90]: The E-I network has an ergodic state where (m 1 = m 2 = 0).Stability analysis to small perturbations of this state reveals an Andronov-Hopf bifurcation towards self-oscillations, when J 11 = 2 and −J 12 J 21 > 1. Upon matching the coefficients of such a linear expansion: we get an approximate mapping into the parameters β, c of the simplest adaptive Ising model: The role of topology One of the most interesting questions about synchronization in neural networks is how general features of the interaction topology affect the collective behavior, i.e., how structure affects function in general terms [91].
From a modeling perspective, most efforts have focused on studying the Kuramoto model (KM) [92], where the individual excitable units are already postulated to be oscillators (for a discussion about this point see [93]).Nevertheless, no exact analytical results for the KM on general networks are available up to now, with an intense debate currently focusing on the nature of the onset of synchronization in strongly heterogeneous topologies [94].
In contrast, we provide here an heuristic argument that the adaptive Ising model directly inherits the wealth of knowledge accumulated about its equilibrium counterpart, in particular, with regard to the features and the location of its critical point(s).The critical point characteristics have been rigorously determined in several geometries, from dimensional lattices to complex networks and small world [95][96][97].
The fact that equilibrium Ising results can be generalized to the adaptive case can be seen directly from the application of the linear response theory, upon considering the Landau expression for the free energy: by construction, the bifurcation point of the dynamical model coincides with the critical point of the underlying equilibrium model.
For instance, for the case of uncorrelated tree-like random graphs, described by a degree distribution P (k) [98], the linear response applied to a Curie-Weiss approximate expression for the free energy [99] leads to the following approximate dynamical equations for the firing rates of nodes with degree k, m k (where k is the mean degree, m = k P (k)m k is the average firing rate and m v = As it can be easily verified by linearizing around the stationary solution m k = h = 0, these equations show that the model has a bifurcation point located at the same position as the equilibrium critical point, i.e. (βJ) c = k 2 k (a more refined calculation [99] based on asymptotically exact Bethe-Peierls approximation gives (βJ )).This simple example shows that the inverse temperature gets renormalized by the branching ratio [98] k 2 k , a topological measure of the density of links, or synaptic connections in our context, that could be considered itself as the key control parameter driving the system in and/or out the synchronized phase.A direct consequence for our case is that if the topology we were considering were scale free, i.e. with an heavy tail for the degree distribution P (k) ∼ k −γ γ < 3, then β c → 0 and the system would always be in the synchronized phase, a feature shared by many collective phenomena in strongly heterogeneous networks [97].In our case, subcritical dynamics is inferred from data and the scale-free topology is not appropriate, but the reasoning here demonstrates clearly how known facts about equilibrium Ising on different topologies can directly translate into the insights of the adaptive Ising model.

APPENDIX C: DETRENDED FLUCTUATIONS ANALYSIS OF THE ALPHA BAND AMPLITUDE ENVELOPE
The DFA [60] consists of the following steps: (i) Given a time series x i (i = 1, ..., N ) calculate the integrated signal I(k) = For a power-law correlated time series, the average r.m.s.fluctuation function F (n) and the box size n are connected by a power-law relation, that is F (n) ∼ n α .The exponent α quantifies the long-range correlation properties of the signal.Values of α < 0.5 indicate the presence of anti-correlations in the time series x i , α = 0.5 absence of correlations (white noise), and α > 0.5 indicates the presence of positive correlations in x i .The DFA was applied to the alpha band (8 − 13 Hz) amplitude envelope.Data were band filtered in the range 8-13 Hz using a FIR filter (second order) designed in Matlab.The scaling exponent α was estimated in the n range corresponding to 2s -60s to avoid spurious correlations induced by the signal filtering [35,52].

APPENDIX D: EXTREME EVENTS, INSTANTANEOUS NETWORK EXCITATION A , NEURONAL AVALANCHES
For each sensor, positive and negative excursions beyond a threshold e were identified.In each excursion beyond the threshold, a single event was identified at the most extreme value (maximum for positive excursions and minimum for negative excursions).Comparison of the signal distribution to the best fit Gaussian indicates that the two distributions start to deviate from one another around ±2.7SD [19].Thus, thresholds smaller than ±2.7SD will lead to the detection of many events related to noise in addition to real events whereas much larger thresholds will miss many of the real events.To avoid noise-related events while preserving most of relevant events, in this study with set the threshold e at ±2.9 standard deviations (SD).The raster of identified events was binned at a number of temporal resolutions , which are multiple of the sampling time T = 1.67 ms.
The network excitation A at a given temporal resolution is defined as the number of events occurring across all sensors in a time bin.An avalanche is defined as a continuous sequence of time bins in which there is at least an event on any sensor, ending with at least a time bin with no events (Fig. 4A).The size of an avalanche, s, is defined as the number of events in the avalanche.For further details see [19,36].

A B
Fig. S3.The reversal time t is defined as the time interval between consecutive zero-crossing events in the ongoing network activity m (Fig. 1).The quantity a0 is the area under the curve between two zero-crossing events.spins, and for MEG data (average over subjects).The model network is parceled in 100 disjoint subsystems, each including 100 spins.In all cases the model is in the resonant regime.The network excitation A is rescaled by the average network excitation A ( = 2 = 2T , where T is the sampling interval, as in Fig. 3).β = 0.99 corresponds to the average β value inferred from MEG data.spins, and for MEG data (average over subjects).The model network is parceled in 100 disjoint subsystems, each including 100 spins.In all cases the model is in the resonant regime.The quiescence duration I is rescaled by the average quiescence duration I ( = 2 = 2T , where T is the sampling interval, as in Fig. 3).β = 0.99 corresponds to the average β value inferred from MEG data.Inset: Probability P0 of quiescence periods as a function of for different β values.We notice that as we move away from the critical point βc = 1, the probability tends to follow an exponential behavior, i.e.P0 ∝ e −a .On the other hand, for β = 0.99 we find P0 ∝ e −a β I , with βI = 0.6304 ± 0.0046, close to the value measured in MEG data (green circles) (βI = 0.5669 ± 0.0117)).

4 DFIG. 1 .
FIG. 1. Adaptive Ising model exhibits coexistance of oscillations and scale-free activity excursions near the critical point.(A) Schematic illustration of the model.Interacting spins si (i = 1, 2, ..., N ) take values +1 (up arrows) or −1 (down arrows) and experience a time-varying external field h(t) that mimics an activity-dependent feedback mechanism.(B) Phase diagram for the mean-field adaptive Ising model.An Andronov-Hopf bifurcation at βc = 1 separates self-sustained oscillations in the total activity m(t) for β > βc (green) from the regime of intermittent oscillations (yellow) for c above c * (β) (solid red line) and an Ornstein-Uhlenbeck process for c below c * (gray).(C) Reversal time t is the time interval between consecutive zero-crossing events in m and a0 is the area under the m(t) curve between two zerocrossing events (inset).Distributions P (a0) are shown in the resonant regime, c > c * , for different values of β.When β ≈ 1, P (a0) is approximately power-law with exponent τ = 1.29 ± 0.01.(D) Distributions P (t) of the reversal times are shown in the resonant regime, c > c * , for different values of β.When β ≈ 1, P (t) is approximately power-law with exponent αt = 1.40 ± 0.01.Inset: Distributions P (a0) and P (t) for the uncoupled model, J = 0, always exhibit exponential instead of power-law behavior (note linear horizontal scale).

FIG. 2 .
FIG. 2. MEG resting state activity of the human brain corresponds to a marginally sub-critical adaptive Ising model.(A) Example trace from a single MEG sensor (top) predominantly contains power in the alpha band (8-13 Hz; bottom, shaded region).Power spectra of MEG signals (bottom; gray = average across 273 MEG sensors for each of the 14 subjects; green = average over sensors and subjects) peak around 10 Hz. (B) Example alpha bandpass filtered MEG signal (top; green) and the simulated total activity m(t) of a model with parameters matched to data (top; violet) show qualitatively very similar behavior.Model parameters (β = 0.9870 and c = 0.1129 for this trace) are inferred by fitting the analytical form of the autocorrelation function C(τ ) (bottom; green line), to autocorrelation estimated from MEG data (bottom; violet dots).(C) Schematic of the area under the curve ae (red) for a given threshold ±e in units of signal SD (top).Distributions P (log ae) of the logarithm of the area under the curve ae, with e = 2.5SD, for MEG data (green curves = average over sensors for each subject) and the model (violet curve = simulation at baseline parameters, see text).Inset: Rescaled distributions of ae collapse to a universal Weibull-like distribution across different threshold values e (Weibull parameters: k = 1.74, λ = 2.58).(D) Central frequency f = ω/2π = (βc − (1 − β) 2 ) 1/2 /8π of the fitted model plotted against fitted β, across all MEG sensors and subjects (color = fraction of total MEG signal power in the alpha band).β values closer to critical βc = 1 are correlated with higher power in the alpha band (R 2 = 0.59; p < 0.001).(E) Root-mean-square fluctuation function F (n) of the DFA for the amplitude envelope of MEG sensor signals in the alpha band (green lines = individual sensors for a single subject).F (n) scales as F (n) ∝ n α for 2 s < n < 60 s (dashed lines), with α > 0.5 for all MEG sensors (0.53 < α < 0.85).(F) Inferred β values correlate with the corresponding DFA exponents α for all MEG sensors and subjects.

FIG. 3 .
FIG. 3. Non-exponential extreme event statistics in MEG resting state activity are reproduced by a marginally subcritical adaptive Ising model.(A) Extreme events on a single sensor are defined as extreme signal excursions (top; red dots) crossing a threshold e = ±nSD (horizontal lines).Resulting raster of extreme events shown across 273 MEG sensors of a single subject across approximately 500 ms recording (bottom).(B) Events are grouped together in temporal bins n = nT in multiples of the sampling interval T (top), to define instantaneous network excitation A , the total number of extreme events across all sensors in a time bin.Representative sequences of network excitation extracted from the raster in the top panel for increasing bin size n (bottom).(C) Rescaled distribution of network excitation, P (A ), for e = 2.9SD and a range of bin sizes n (different plot symbols) in MEG data (green symbols; average over subjects) and in the model simulated at baseline parameters with K = 100 subsystems of n sub = 100 neurons each (violet symbols).Distributions for different collapse onto a single non-exponential master curve for both data and model.Corresponding distribution in phase-scrambled MEG signals shows an exponential behavior, with absence of high excitation events (brown = surrogate data).Inset: Rescaled P (A ) (green = average over subjects; violet = average over model simulations) and respective standard deviation (colored areas) shown for bin size = 2T.(D) Rescaled distributions of quiescence durations, P (I ) collapse onto a single master curve for different .Plotting conventions and model simulation details are the same as in (C).Top inset: Rescaled P (I ) (green = average of subjects; violet = average over model simulations) and respective standard deviation (colored area) shown for bin size = 2T.Bottom inset: Probability P0 of finding a quiescent time bin scales approximately as P0 = exp −a β I with bin size ; βI = 0.582 ± 0.013 and βI = 0.610 ± 0.012 for data and model, respectively; βI = 0.996 ± 0.001 for surrogate data.
FIG. 4. Scale-free neuronal avalanches in MEG resting state activity are reproduced by a marginally subcritical adaptive Ising model.(A) Schematic representation of a neuronal avalanche.Avalanche size s is the sum of network excitations A over time bins belonging to the avalanche; its duration d is the number of bins times their duration, .(B) Distribution of avalanche sizes, P (s), for MEG data (green circles with error bars = average over subjects ± standard deviation) and the model simulated at baseline parameters with K = 100 subsystems of n sub = 100 neurons each (violet squares with error bars = average over model simulations ± standard deviation).Both distributions are estimated using a threshold e = 2.9SD and bin size 4 = 4T .The brown curve is the distribution P (s) obtained from surrogate data (Appendix E) with the same threshold and bin size.Inset: Average avalanche size scales with its duration as s ∼ d ζ (different plot symbols = different as in Fig. 3; green = MEG data; violet = model simulation; model simulation curves are vertically shifted for clarity), so that the exponent ζ remains independent of the time bin size .ζ = 1.28 ± 0.01 for MEG data (dashed line) and ζ = 1.58 ± 0.03 for model simulation (thick line).
k i=1 (x(i) − x ), where x is the mean of x i ; (ii) Divide the integrated signal I(k) into boxes of equal length n and, in each box, fit I(k) with a first order polynomial I n (k), which represents the trend in that box; (iii) For each n, detrend I(k) by subtracting the local trend, I n (k), in each box and calculate the root-mean-square (r.m.s.) fluctuation F (n) = N k=1 [I(k) − I n (k)] 2 /N ; (iv) Repeat this calculation over a range of box lengths n and obtain a functional relation between F (n) and n.

Fig. S1 .Fig. S2 .
Fig. S1.Autocorrelation C and ongoing network activity m for β = 0.5 and different c values.Far from the critical point, the presence of a strong adaptive feedback may also produce short -C rapidly decays to zero -intermittent oscillation bursts.(A) Autocorrelation for different c values.(B) m for c = 0.5.(C) m for c = 2. (D) m for c = 10.tc is the inferred autocorrelation time from exponential fit.

Fig. S4 . 99 Fig. S6 .Fig. S7 .Fig. S8 .
Fig.S3.The reversal time t is defined as the time interval between consecutive zero-crossing events in the ongoing network activity m (Fig.1).The quantity a0 is the area under the curve between two zero-crossing events.(A) Distribution P (a0) of the quantity a0 for the model at the critical point β = 1 for the different strengths c of the adaptive feedback.(B) Distribution P (t) of the reversal time for β = 1 and different values of the parameter c.

Fig. S9 .Fig. S10 .
Fig. S9.Distributions of the activity per bin A in model simulations with N = 10 4 and N = 9 • 10 4 spins (β = 0.99, c = 0.01).The model network is parceled into subsystems of different size n sub .The distribution P (A ) is independent of the subsystem size n sub .The network excitation A is rescaled by the average network excitation A , with = 2T , where T is the sampling interval as in Fig. 3.

N/n sub = 200 AFig. S11 . 99 Fig. S13 .
Fig. S11.Distribution P (A ) (A), P (I ) (B), and P (s) (C) in model simulations with N = 9 • 10 4 spins (β = 0.99, c = 0.01) and different numbers of subsystems K = N/n sub .(A) The distribution of network excitation P (A ) ( = 2 = 2T ) weakly depends on the number of subsystems K = N/n sub .Inset: Distributions P (A ) for different values of the threshold e used to detect extreme events.The number of subsystems is fixed to N/n sub = 200.(B) The distribution of quiescence durations P (I ) ( = 2 = 2T ) depends on the number of subsystems K = N/n sub , particularly on the tail (black triangles up).For N/n sub = 200 a good agreement between data and model simulations is recovered when the threshold e is increased from 2.9SD (the value used for N/n sub = 100) to 3.3SD (orange triangles up).Inset: Probability P0 of finding a quiescent time bin scales approximately as P0 = exp −a β I with bin size ; βI 0.6 depends on the number of subsystems N/n sub , and slightly increases from ≈ 0.6 for N/n sub = 100 (violet squares) to ≈ 0.7 for N/n sub = 200 (black triangles).(C) The distribution of avalanche sizes P (s) weakly depends on the number K = N/n sub .Inset: Distributions P (s) for different values of the threshold e used to detect extreme events.The number of subsystems is fixed to N/n sub = 200.Increasing e has the effect of decreasing (increasing) the probability of large (small) avalanches.All P (s) distributions are calculated with = 4 = 4T .

Fig. S14 .
Fig. S14.Distributions of quiescence durations I for different values of β in model simulations with N = 10 4spins, and for MEG data (average over subjects).The model network is parceled in 100 disjoint subsystems, each including 100 spins.In all cases the model is in the resonant regime.The quiescence duration I is rescaled by the average quiescence duration I ( = 2 = 2T , where T is the sampling interval, as in Fig.3).β = 0.99 corresponds to the average β value inferred from MEG data.Inset: Probability P0 of quiescence periods as a function of for different β values.We notice that as we move away from the critical point βc = 1, the probability tends to follow an exponential behavior, i.e.P0 ∝ e −a .On the other hand, for β = 0.99 we find P0 ∝ e −a β I , with βI = 0.6304 ± 0.0046, close to the value measured in MEG data (green circles) (βI = 0.5669 ± 0.0117)).

70 Fig. S15 .
Fig. S15.Distribution of avalanche sizes, P (s), for MEG data (green curve = average over subjects) and the model simulated at different β values in the resonant regime, i.e. c > c * .Distributions are estimated using a threshold e = 2.9SD and bin size 4 = 4T .Already for β = 0.95, a value slightly smaller than the baseline value 0.99, avalanche sizes from model simulations tend to follow an exponential distribution that is far from reproducing avalanche size distributions from MEG data.
Fig. S16.Distribution P (s) of avalanche sizes in model simulations with N = 10 4 and N = 9 • 10 4 spins (β = 0.99, c = 0.01).The model network is parceled in subsystems of different size n sub .The distribution P (s) is independent of the subsystem size n sub .