## Abstract

Criticality has gained widespread interest in neuroscience as an attractive framework for understanding the character and functional implications of variability in brain activity. The metastability of critical systems maximizes their dynamic range, storage capacity, and computational power. Power-law scaling—a hallmark of criticality—has been observed on different levels, e.g., in the distribution of neuronal avalanches *in vitro* and *in vivo*, but also in the decay of temporal correlations in behavioral performance and ongoing oscillations in humans. An unresolved issue is whether power-law scaling on different organizational levels in the brain—and possibly in other hierarchically organized systems—can be related. Here, we show that critical-state dynamics of avalanches and oscillations jointly emerge in a neuronal network model when excitation and inhibition is balanced. The oscillatory activity of the model was qualitatively similar to what is typically observed in recordings of human resting-state MEG. We propose that homeostatic plasticity mechanisms tune this balance in healthy brain networks, and that it is essential for critical behavior on multiple levels of neuronal organization with ensuing functional benefits. Based on our network model, we introduce a concept of multi-level criticality in which power-law scaling can emerge on multiple time scales in oscillating networks.

## Introduction

Synchronous brain activity is thought to be crucial for neural integration, cognition, and behavior (Fries, 2005; Buzsáki, 2006). The multi-scale properties of synchronous cell assemblies, however, remain poorly understood. While probing activity at different scales, several investigators have begun to consider self-organized criticality (SOC) as an overriding neuronal organizing principle (Linkenkaer-Hansen et al., 2001; Beggs and Plenz, 2003; Haldeman and Beggs, 2005; Kinouchi and Copelli, 2006; Levina et al., 2007; Plenz and Thiagarajan, 2007; Gireesh and Plenz, 2008; Poil et al., 2008; Petermann et al., 2009; Chialvo, 2010; He et al., 2010; Millman et al., 2010; Rubinov et al., 2011). The SOC theory holds that slowly driven, interaction-dominated threshold systems will be attracted to a critical point, with the system balanced between order and disorder (Bak et al., 1987; Bak, 1996; Jensen, 1998).

This critical state is characterized by scale-free probability distributions. In cortical slices, neuronal avalanches of local field potential activity are governed by a power-law (i.e., scale-free) distribution of event sizes with an exponent of −1.5, as expected for a critical branching process (Beggs and Plenz, 2003). Similar results have also been found *in vivo* (Petermann et al., 2009). This finding inspired the development of computational models that were capable of producing power-law-distributed neuronal avalanches through a process of self-organization, thereby providing theoretical support for SOC in neuronal systems (Levina et al., 2007; Millman et al., 2010; Rubinov et al., 2011). Models have shown that networks in a critical state display differing responses to the widest input range (Kinouchi and Copelli, 2006), with the largest repertoire of metastable patterns (Haldeman and Beggs, 2005), and maximal information and transmission capacity (Beggs, 2008; Shew et al., 2011). Scale-free dynamics have been observed on many levels of neuronal organization. Power-law scaling is, e.g., found in the decay of temporal correlations in behavioral performance (Gilden, 2001) and in the phase-locking intervals between brain regions (Gong et al., 2003; Kitzbichler et al., 2009).

Interestingly, ongoing oscillations measured with EEG or MEG exhibit power-law-scaled temporal (auto)correlations in their amplitude modulation, also known as long-range temporal correlations (Linkenkaer-Hansen et al., 2001; Monto et al., 2007). Given that neuronal avalanches and oscillations both depend on balanced excitation and inhibition (Beggs and Plenz, 2003; Atallah and Scanziani, 2009; Shew et al., 2011), it is plausible that their scale-free dynamics are related. To investigate this, we modeled a generic network of excitatory and inhibitory neurons in which the local connectivities could be varied systematically (Fig. 1). We found that scale-free dynamics of avalanches and oscillations jointly emerge from balanced excitatory and inhibitory connectivity. We introduce the concept of multi-level criticality as a state where scale-free behavior emerges on different levels of network dynamics: the short-time-scale spreading of activity, with an upper bound at the characteristic time scale of the dominant oscillation, and the long-time-scale modulation of the oscillatory amplitude.

## Materials and Methods

##### Network organization.

We modeled networks of excitatory (75%) and inhibitory (25%) neurons (*n* = 2500) arranged in an open grid, with local functional connectivity (Fig. 1*A*,*B*). To define the local connectivity, each neuron was given a square area (width = 7 neurons; Fig. 1*A*) where they could connect to 48 other neurons. The excitatory connectivity was set between 10–70% and the inhibitory connectivity was set between 30–90% of the total number of neurons within this local range. The neurons were set on a 2D open grid. Border neurons had fewer connections, because these neurons had a lower number of neurons in their local range. In neocortical tissue, local functional connectivity between neurons decays with distance; however, only a few empirical reports have assessed this quantitatively (Holmgren et al., 2003; Sporns and Zwi, 2004). We therefore—given the small range of connectivity in our model—implemented the connectivity with an exponential decay. More specifically, the probability, *P*, of a connection at a distance *r* was given by *P*(*r*) = *Ce*^{−r}, where *C* is a constant determining the connectivity probability.

We created five networks for each combination of excitatory and inhibitory connectivity, and each of these networks were simulated for 1000 s. The weights (*W _{ji}*) of these connections were dependent on the type (excitatory or inhibitory) of the presynaptic (

*i*) and postsynaptic (

*j*) neuron (Fig. 1

*B*), but were otherwise fixed to limit the number of free model parameters (for specific weights values, see Neuron model, below). The models implemented randomness at the level of assigning whether a neuron was excitatory or inhibitory, and in deciding whether neurons would connect to each other.

##### Neuron model.

Neurons were modeled using a synaptic model integrating received spikes, and a probabilistic spiking model (Fig. 1*C*). Each time step (*dt*) of 1 ms starts with each neuron, *i*, updating *I _{i}* with received input from any of the set of connected neurons (

*J*), together with an exponential synaptic decay, which can either be excitatory or inhibitory.

*S*is a binary spiking vector of the neurons spiking in previous time step, and

*W*= 0.011,

_{IE}*W*= 0.02,

_{EE}*W*= −2, and

_{EI}*W*= −2. τ

_{II}_{I}= 9 ms and

*I*

_{0}= 0 for both excitatory and inhibitory inputs. The probability of spiking,

*P*is then updated with this input, together with an exponential decay: with τ

_{S}_{P}(excitatory) = 6 ms, τ

_{P}(inhibitory) = 12 ms,

*P*

_{0}is the background spiking probability with

*P*

_{0}(excitatory) = 0.000001 [1/ms],

*P*

_{0}(inhibitory) = 0 [1/ms].

We determine whether the neuron spikes with the probability *P _{S}*, and update the spiking vector for the next time step. If a neuron spikes, the probability is reset to the reset value

*P*(excitatory) = −2 [1/ms] and

_{r}*P*(inhibitory) = −20 [1/ms], and the binary spiking vector

_{r}*S*is updated. In the next time step, all neurons that it connects to will have their input updated according to Equation 1.

##### Human MEG.

Ongoing brain activity was measured with MEG in sessions of 20 min from eight normal subjects (aged 20–30 years, 2 females). The subjects were seated in a magnetically shielded room and instructed to relax and sit still with eyes closed during the recording. The data were recorded with 204 planar gradiometers, sampled at 900 Hz, and decimated off-line to 300 Hz with a passband of 0.1–100 Hz (sixth-order Butterworth digital filters). The data have been used in Linkenkaer-Hansen et al. (2004) and the study was approved by the Ethics Committee of the Department of Radiology of the Helsinki University Central Hospital.

##### Detrended fluctuation analysis of long-range temporal correlations.

The detrended fluctuation analysis (DFA) was used to analyze the scale-free decay of temporal (auto)correlations, also known as long-range temporal correlations (LRTC). The DFA was introduced as a method to quantify correlations in complex data with less strict assumptions about the stationarity of the signal than the classical autocorrelation function or power spectral density (Linkenkaer-Hansen et al., 2001). An additional advantage of DFA is the greater accuracy in the estimates of correlations, which facilitates a reliable analysis of LRTC up to time scales of at least 10% of the duration of the signal (Chen et al., 2002; Gao et al., 2006). DFA exponents in the interval of 0.5 to 1.0 indicate scale-free temporal correlations (autocorrelations), whereas an exponent of 0.5 characterizes an uncorrelated signal. The main steps from the broadband signal to the quantification of LRTC using DFA have been explained in detail previously (Linkenkaer-Hansen et al., 2001). In brief, the DFA measures the power-law scaling of the root-mean-square fluctuation of the integrated and linearly detrended signals, *F*(*t*), as a function of time window size, *t* (with an overlap of 50% between windows). The DFA exponent is the slope of the fluctuation function *F*(*t*), and can be related to the power-law scaling exponent of the autocorrelation function.

##### Peak frequency analysis.

To estimate the peak frequency, we fitted a Gaussian to the power spectrum estimate (Welch method) in the expected peak frequency range.

##### Quantification of power-law scaling using the κ index.

We visually inspected the avalanche distributions and applied both least-square fitting, as well as maximum likelihood method fitting (Clauset et al., 2009) to those showing log–log linear scaling. This preliminary analysis pointed to power-law exponents of −1.5, in agreement with previous literature on neuronal avalanches (Beggs and Plenz, 2003; Mazzoni et al., 2007; Gireesh et al., 2008; Pasquale et al., 2009; Petermann et al., 2009; Shew et al., 2011; Yang et al., 2012). However, our model clearly did not produce log–log linear scaling for all parameter combinations and, therefore, we could not use maximum likelihood fitting and power-law exponents to identify subcritical or supercritical dynamics. Instead, we quantified the similarity between the distribution of our data and a power-law by calculating the average difference of the cumulative distribution of a power-law function with an exponent of −1.5 and that of our experimental data at 10 equally spaced points on a logarithmic axis and adding one (Fig. 2*D*). This gives the κ index (Shew et al., 2009). A subcritical distribution is characterized by κ < 1, and supercritical distribution by κ > 1, whereas κ ≈ 1 indicates a critical network. Using the κ index had the added advantage of making our results comparable to several recent studies using this index (Shew et al., 2009, 2011).

## Results

The model displayed many forms of activity propagation depending on the excitatory/inhibitory connectivity ratio (*E*/*I* balance). This is illustrated in Figure 2*A*, where the temporal dynamics of activity propagation can be inferred from the color scale indicating the time since a given neuron spiked. For the networks with a low *E*/*I* balance, the activity was characterized by localized wave-like spreading (Fig. 2*A*, left). The propagation stops due to the strong local inhibition. With increasing *E*/*I* balance, the waves were able to spread further until they reached across the entire network. In these networks, sustained activity patterns were possible (Fig. 2*A*, middle), with complex patterns such as spiral waves repeating themselves several times. These patterns would eventually dissipate with the network either going through a period of low activity similar to the low *E*/*I* balance case, or the network would become highly active and lose any spatially contiguous patterns. This highly active state is characteristic of the networks with a high *E*/*I* balance (Fig. 2*A*, right).

To quantify the propagation of activity in the network, we investigated whether it could be reconciled with the concept of neuronal avalanches reported for local field potential recordings *in vitro* (Beggs and Plenz, 2003) and *in vivo* (Petermann et al., 2009). To this end, we defined a neuronal avalanche using a modified version of the method developed by Beggs and Plenz (2003). As some of the networks analyzed were continuously active, we altered their method to define the start and end of an avalanche by a time step with spatially integrated network activity below the threshold value set to 50% of the median spike activity (Fig. 2*B*). The size of an avalanche is the total number of spikes during the avalanche. We characterized the size distributions (Fig. 2*C*) using the statistical measure κ (Fig. 2*D*), which was developed by Shew et al. (2009) (see Material and Methods, above). In agreement with visual inspection of the distributions, κ indicated that the model could produce subcritical, critical, or supercritical dynamics depending on the connectivity parameters (Fig. 2*C*,*E*). The subcritical regime (κ < 1) had exponentially decaying avalanche sizes, the supercritical regime (κ > 1) had a characteristic scale indicated by the bump at large avalanche sizes, and power-law scaling extended up to the system size in the transition regime (κ ≈ 1).

On longer timescales (>120 ms), the progression of these avalanches gave the spatially integrated network activity a distinct oscillatory character (Fig. 3*A*), which was qualitatively similar to what is typically observed in recordings of human resting-state MEG (Fig. 3*B*,*C*) (Linkenkaer-Hansen et al., 2004; Poil et al., 2008). Most network configurations produced 8–16 Hz oscillations (Fig. 3*D*,*E*); however, a clear transition line from low to high power is visible in Figure 3*F*. This transition line marks the shift in network activity from sparse activity with low synchrony to strongly oscillating activity, while activity in between these regimes is characterized by waxing and waning oscillations.

To quantitatively assess scaling behavior on the level of oscillations, we analyzed the waxing and waning of the oscillations (Fig. 3*B*) on the long time scales of 2–10 s using detrended fluctuation analysis (DFA) (Fig. 3*G*) (Linkenkaer-Hansen et al., 2001). DFA characterizes the power-law scaling of the temporal autocorrelations of a signal, so-called long-range temporal correlations (see Materials and Methods, above). If the oscillations wax and wane at random (e.g., as bandpass filtered white noise), the DFA exponent will be 0.5, whereas a DFA exponent between 0.5 and 1.0 imply that oscillations have a scale-free amplitude modulation (Linkenkaer-Hansen et al., 2001). Large-scale studies have reported DFA exponents of alpha-band oscillations in the range of 0.55–1.05 for adults (Smit et al., 2011); however, most often exponents fall in the range of 0.6–0.9, which is in agreement with the empirical data included in this study (α = 0.82 ± 0.09 SD, *n* = 9, average across the four parietal electrodes with the highest amplitude). In our model, DFA revealed a critical regime of strong temporal correlations (α = 1.0), whereas a relative increase in either inhibitory or excitatory connectivity led to reduced temporal correlations (α = 0.6), similar to those of filtered white noise (α = 0.5; Fig. 3*H*). Of note, variability of DFA exponents between networks with the same connectivity ratio was large and highest in the critical region (Fig. 3*I*), suggesting that individual differences in long-range temporal correlations in humans may be explained by subtle changes in cortical networks operating near criticality. Interestingly, the binned averages of DFA exponents peak when the network produces critical avalanches, with the correlations decreasing for the more subcritical or supercritical avalanches (Fig. 4). Together, our data indicate that nontrivial scaling behavior in the modulation of oscillatory activity can emerge on much longer time scales than those implemented in the model, and that neuronal avalanches may be viewed as the building blocks of neuronal oscillations with scale-free amplitude modulation.

## Discussion

We have investigated the dynamics of oscillatory neuronal networks and discovered that the networks produced scale-free dynamics of neuronal avalanches (<100 ms) and oscillations (>2 s) through a common mechanism of balanced excitatory and inhibitory connectivities. We propose a concept of multi-level criticality to denote the phenomenon of critical behavior on multiple levels. Multi-level criticality may define a new class of dynamical systems, because it allows for criticality to emerge jointly on multiple levels separated by a characteristic scale, which is traditionally considered contradictory for critical dynamics. In the following, we comment on the model parameters and how our observations relate to previous reports on neuronal avalanches, long-range temporal correlations in oscillations, and the predictions that can be derived.

### Model choices and mechanisms of spatial and temporal correlations

Interaction between excitatory and inhibitory neurons is important for the generation of oscillations (Buzsáki and Draguhn, 2004) and neuronal avalanches (Beggs and Plenz, 2003). Therefore, we developed a model in which variation in excitatory and inhibitory connectivity could be systematically and independently varied. Other parameters were kept fixed to limit the number of parameter combinations that required investigation (see Materials and Methods, above).

Local connectivity was implemented, because it is observed empirically (Holmgren et al., 2003; Sporns and Zwi, 2004) and supports the emergence of oscillations and propagating patterns, both of which were important for the present study. For example, we would not have observed the spreading of activity in Figure 2*A* with a random network topology. The restricted spreading of activity is essential for build up of temporal autocorrelations and criticality, because it allows blobs of nonexcitable and excitable neurons to form. These blobs act as memories of past activity, shaping the path of future activity, analogous to a river's flow being shaped by the past flow of water (Bak, 1996; Chialvo and Bak, 1999). Had it not been for the stochastic spiking implemented in our model, the activity of the network would follow more stereotyped paths, such as persistent oscillations. A recent study of model networks showed that a hierarchical topology can increase the range of parameters that give rise to critical dynamics—and thereby also the robustness of a critical state—because the modular topology limits activity spreading to local parts of the network (Wang and Zhou, 2012). Similarly, our localized topology is likely to support the emergence of oscillations with long-range temporal correlations by allowing the build up of multiple spatially contiguous patterns that repeat over many oscillation cycles. Thus, the specific network topology and the random spiking in the present model neuronal networks are likely to have influenced the emergence of criticality, which could be the subject of future investigations.

### Linking model behavior to empirical observations

To our knowledge, there are no empirical reports investigating neuronal avalanches and long-range temporal correlations from the same dataset. However, that characteristic scales of oscillations are compatible with the scale-free spreading of activity in neuronal avalanches is supported by *in vivo* and *in vitro* data (Gireesh et al., 2008). Similar to our model, this study reported that propagation of local field potential activity during cycles of nested theta (4–15 Hz), beta (15–30 Hz), and gamma (30–100 Hz) oscillations showed power-law scaling.

Our model produces a large variation of DFA exponents in the transition region to criticality (Fig. 3*H*,*I*), and also a higher variability of DFA and κ index in the critical region (Figs. 2*F*, 3*I*). Several EEG/MEG studies have reported DFA scaling exponents in the alpha-frequency band (8–13 Hz) to be ∼0.75 (Linkenkaer-Hansen et al., 2001; Nikulin and Brismar, 2004). Underlying this mean value is a large interindividual variability from ∼0.55 to 1.05, as shown in a recent large-scale study involving >1400 subjects (Smit et al., 2011). Our model suggests that subtle differences in the balance between excitation and inhibition could explain the empirically observed variation in DFA exponents. Still, the interesting question of whether empirically observed long-range temporal correlations emerge from a slightly subcritical or supercritical network remains. Here, we notice that a κ index between 1.02 and 1.14 has been reported for neuronal avalanches (Shew et al., 2011; Yang et al., 2012), which in our model would correspond to a DFA exponent in the range of 0.9 to 0.65 (Fig. 4*B*) and points to oscillation-generating networks being slightly toward the supercritical regime.

The higher variability of DFA and κ index in the critical region (Figs. 2*F*, 3*I*) also correlates well with the findings of higher entropy and phase variability in empirical studies of neuronal avalanches (Shew et al., 2011; Yang et al., 2012) and studies showing maximal dynamic range and number of metastable states in this region (Haldeman and Beggs, 2005; Kinouchi and Copelli, 2006). This variability reflects the inherent dynamics of a critical system, which lives on the tension zone between two distinct states (Bak, 1996; Christensen and Moloney, 2005). A potential connection between long-range temporal correlations and variability in phase synchrony could be that a network with maximal variability in phase synchrony can better maintain scale-free dynamics with a wide dynamic range.

In view of the subtle changes in parameters of the model that cause changes in the DFA and the considerable variation in network topologies underlying oscillations in different frequency bands and cortical regions (Miller et al., 2008), it may be premature to propose an optimal DFA exponent. For example, preclinical studies have shown a decrease in the DFA of alpha oscillations in parietal regions of Alzheimer patients (Montez et al., 2009), an increase in the beta oscillations in epileptogenic zones (Monto et al., 2007), and schizophrenia patients show strongly attenuated LRTC in both alpha and beta oscillations in widespread cortical regions (Nikulin et al., 2012). To better understand the functional implication of LRTC in oscillations, future work should address, e.g., whether information transfer is optimal in a multilevel critical network, as it has been shown in a system producing neuronal avalanches (Haldeman and Beggs, 2005; Kinouchi and Copelli, 2006).

### Predictions derived from the model

The main prediction derived from our model is that avalanche dynamics and long-range temporal correlations are related, which to our knowledge has never been addressed empirically. Specifically, the model predicts that the dominant oscillation will have the strongest LRTC when the activity is in the form of neuronal avalanches, as defined by a κ index of 1. Testing this hypothesis would presumably require invasive sampling of local field potentials *in vivo* using 2D or 3D electrode arrays, and analyzing the data both for neuronal avalanches and long-range temporal correlations (Sirota et al., 2008; Goutagny et al., 2009; Petermann et al., 2009).

Pharmacological manipulation of GABAergic neurotransmission has been shown to influence the strength of LRTC in the beta-frequency range in humans (Monto et al., 2007) and *in vitro* (Poil et al., 2011), and to abolish the power-law scaling of neuronal avalanches (Beggs and Plenz, 2003; Gireesh et al., 2008; Yang et al., 2012), which also follows from our model. Ideally, future studies will perturb the dynamics of cortical networks using different tasks or pharmacological manipulation to validate the correspondence between loss of criticality at the level of neuronal avalanches (seen in a change of the κ index) and at the level of oscillations (seen in lower DFA exponents) within one experiment. In view of the nonlinear relation between the κ index and DFA exponents, it is not possible to rigorously predict the covariation of these indices. In fact, investigating this issue is an important step in deciding whether neuronal networks generating network oscillations with long-range temporal correlations are closest to the subcritical or the supercritical regime.

The present model may pave the way for an integrated understanding of neuronal avalanches and oscillations in information processing and their alteration with disease (Cassenaer and Laurent, 2007; Montez et al., 2009; Lisman, 2010; Shew et al., 2011; Carandini, 2012). As such, future studies should investigate a greater regime of model parameters, e.g., larger networks, different connectivities, or dynamical synapses (Levina et al., 2007; Rubinov et al., 2011). Several mechanisms have, for instance, been shown to support homeostatic plasticity by which synaptic strengths are dynamically adjusted to promote network stability (Turrigiano and Nelson, 2004; Stewart and Plenz, 2006) and could, thus, also be essential for critical-state dynamics. Of note, these mechanisms are dysfunctional in a variety of neuropsychiatric disorders—including autism and schizophrenia—causing changes to the excitatory/inhibitory balance of neuronal microcircuitry and ensuing behavioral impairments (Yizhar et al., 2011).

## Footnotes

This work was supported by Netherlands Organization for Scientific Research (NWO) Physical Sciences Grant 612.001.123 to K.L.-H., a talent grant to K.L.-H. (CvB, VU University Amsterdam), a Young Talent grant to R.H. from Neuroscience Campus Amsterdam, and a Top Talent grant to S.-S.P. (NWO). We thank Dante Chialvo for discussions and two anonymous reviewers for constructive comments.

The authors declare no financial conflicts of interest.

- Correspondence should be addressed to Dr. Klaus Linkenkaer-Hansen, Department of Integrative Neurophysiology, Center for Neurogenomics and Cognitive Research, VU University Amsterdam, De Boelelaan 1085, 1081 HV Amsterdam, The Netherlands. klaus.linkenkaer{at}cncr.vu.nl