M ar 2 01 8 Two types of criticality in the brain

David Dahmen, Sonja Grün, 2 Markus Diesmann, 3, 4 and Moritz Helias 4 Institute of Neuroscience and Medicine (INM-6) and Institute for Advanced Simulation (IAS-6) and JARA Institute Brain Structure-Function Relationships (INM-10), Jülich Research Centre, Jülich, Germany Theoretical Systems Neurobiology, RWTH Aachen University, Aachen, Germany Department of Psychiatry, Psychotherapy and Psychosomatics, Medical Faculty, RWTH Aachen University, Aachen, Germany Department of Physics, Faculty 1, RWTH Aachen University, Aachen, Germany (Dated: July 22, 2018)

The brain is a dynamical system with potentially different regimes of operation. Network models and experiments suggest optimal computational performance at critical points which mark the transition between two dynamical regimes [1]. At critical points, systems exhibit universal behavior, characterized by strong concerted fluctuations between all constituents leading to effective long-range interactions despite short-range couplings.
A particular type of criticality occurs in structurally balanced neuronal networks with excitatory and inhibitory feedback on par, leading to neuronal avalanches [2]. Signatures of avalanches have been observed in indirect measures of neuronal activity in macaque motor cortex [3]. Such evidence, however, is harder to obtain from direct measures of parallel neuronal activity; for extended periods large transients of population activity and long-range temporal correlations are absent (Fig. 1). The experimental data rather show weak and fast fluctuations of the population activity, suggesting an excess of inhibitory feedback and dynamically balanced activity [4]. Average correlations are low as a consequence [5,6]. Experimental evidence for the operation of cortical networks in such balanced states is overwhelming [7,8,10,11]. Such states are also consistent [12] with the observation of wave-like activity [13].
Another type of criticality, devoid of avalanches, has been investigated in computational studies [15]: edge-of-chaos criticality. With increasing heterogeneity in network connections, the dynamics changes from regular to chaotic, as characterized by the largest Lyapunov exponent. Such states moreover possess rich dynamics: a collection of coexisting network modes leads to topologically complex responses [16]. The onset of chaotic dynamics coincides with the breakdown of linear stability of the deterministic dynamics [17]. But a direct measurement of the Lyapunov exponent is hampered by noise and the incomplete knowledge of the dynamical state; as is the distance from the breakdown of linear stability.
We therefore need a measure that indirectly assesses either of these characteristics. Correlations in the neuronal activity ( Fig. 1) may be informative, because they arise from interactions between constituents, which cause the collective behavior. Pairwise covariances c(W ) as shown in Fig. 1, to leading order follow a unique and model invariant simple law [18][19][20] relating them to the effective connectivity matrix W of the linearized dynamics. The covariance matrix D measures fluctuations of the single-unit activity, determined by firing rates as well as by shared and correlated external inputs.
The main contribution to covariances c(W ) arises from interactions within the local network [21]. Therefore, the experimentally observed variability across neurons (Fig. 1D) is predominantly caused by the structural variability contained in W . A potential bias in the estimation of covariances due to a finite amount of data [8] can be accounted for [9]. The linearized network dynamics corresponding to Eq. (1) becomes unstable if one eigenvalue of W has a positive real part (Re(λ) ≥ 1). Eigenvalues are, in principle, determined by all connections in the network. But their determination from experimentally observed covariances by Eq. (1) is hindered by severe subsampling; even with present massively parallel recording techniques at most hundreds of neurons can be measured simultaneously from the same local network [14]. We therefore need to find a characteristics of correlations that is informative about the eigenvalues, insensitive to the details of the network, and that can be estimated from a few hundred neurons at most.
The distribution of covariances is such a candidate measure. We first study how it is affected by the strength of connections (Fig. 2). With increasing connection strength, activity successively propagates over multiple synapses, leading to indirectly mediated interactions via a growing number of parallel paths ( Fig. 2A) [19,20]. The eigenvalues of the connectivity matrix approach the critical line where the dynamics loses linear stability (Fig. 2B); distributions of covariances (Eq. 1) become monomodal and broader but stay centered approximately around zero. These results expose a unique hallmark of dynamically balanced networks that operate close to linear instability: widely distributed covariances with a small mean as observed in motor cortex (Fig. 1D).
Although local cortical networks show non-random, cell-type specific, and distance-dependent connectivity, the simple model of a homogeneous random network studied in Fig. 2 is sufficient to explain gross features of the experimental data. To expose the fundamental mechanism explaining the dispersion of correlations, we change the perspective: instead of treating a single realization of the connectivity and calculating the variance of c ij across pairs, we fix one pair and average over realizations of networks; a well-established technique for disordered physical systems. The disordered connectivity appears as an inverse matrix in Eq. (1), which technically complicates the analysis: no results from random matrix theory apply for this particular problem. Instead we construct a moment generating function [22] for the linearized network dynamics [9], which allows the use of spin-glass techniques [23,24] combined with approximations for large-N field theories [25]. As a result, the disorder contained in the ∼ N 2 entries of the connectivity formally reduces to only two fluctuating auxiliary variables that provide input to a fully symmetric all-to-all connected network. This high symmetry and the drastic reduction of dimensionality yield a mean-field theory beyond self-averaging [9]. The mean and standard deviation of variances (i = j) and covariances (i = j), to leading-order in the network size N , then follow as and synaptic weight w where µ ij ∼ O(1/ √ N ) is the mean and σ 2 ij = λ max /N ∼ O(1/N ) the variance of connection weights in W . The latter determines the radius λ max of the cloud of eigenvalues (Fig. 2B) and the renormalized covariance D λ = D/(1 − λ 2 max ), which accounts for the structural variability of connections (Fig. 3A,B). Eq. (2) explains the low mean c ij ∼ O(1/N ) of covariances as a consequence of the dependence on µ < 0: activity is decorrelated via excess inhibition [5,6]; Eq.
The theory in principle determines the largest eigenvalue from measured covariances. But there are two complications: D λ is not known and cannot be measured, and it is unclear how robust the result is with regard to more realistic connectivity. Both problems are solved by considering the normalized width of the distribution of covariances ∆ = δc ij /c ii : This measure is predominantly determined by the network size and the most unstable eigenvalue λ max . The prediction of ∆ is sufficient even for such network topologies (Fig. 3C) as excitatory-inhibitory networks (Fig. 3F) and distance-dependent connection probabilities (Fig. 3G). The latter networks also qualitatively explain the shape of the experimentally observed covariance distribution. The applicability of the theory goes beyond linear network dynamics; it predicts ∆ even in spiking networks (Fig. 3H). Therefore, ∆ can be used to infer the operational regime of the cortical network (Fig. 1). The distance to linear instability is determined by λ max which, to leading order in N , is given by inversion of Eqs. (2) and (3) as (4) Fig. 3D suggests that given the experimentally measured distribution of covariances (Fig. 1D) and biologically plausible N above 10 4 neurons, the network needs to operate close to instability to explain the data. The weak dependence of λ max on N and ∆ in the biologically relevant ranges further secures the robustness of the result towards variations of these experimental estimates. In summary, we combine methods from field theory and disordered systems to predict the dispersion of neuronal covariances in relation to the operational regime of the network. The theory goes beyond the mean-field approximation of population-averaged dynamics and, applied to massively parallel spike recordings from macaque motor cortex, unravels its dynamics close to criticality.
The dynamics in this regime is in contrast to strongly fluctuating population activity observed in avalanches [26]. These cause positive average covariances ( Fig 4B) and a slowly decaying autocorrelation function of the population activity (Fig. 4C). A close structural balance between excitation and inhibition (µ ≈ 1/N ) explains such behavior, caused by a single, nearly unstable eigenvalue that amplifies fluctuations of the population activity [27] (Fig 4A, red dot). All other modes have low amplitude and the same exponentially decaying, fast dynamics ( Fig. 4C inset). The network is hence effectively one-dimensional [28].
In general, each eigenvalue corresponds to a specific linear combination of neuron activities, a mode ( Fig. 4A, D). The deterministic network structure (structural balance, populations, hubs, etc) determines only a small subset of eigenvalues, whereas heterogeneity of connections across neuron pairs shapes the bulk. If sufficiently large, this variability causes a multitude of modes close to instability (Fig. 4D). The broad distribution of covariances in the cortical data (Fig. 1D, Fig. 4E) implies such a critical state. The almost vanishing mean covariances and the weakly fluctuating and quickly decaying population activity (Fig. 1F, Fig. 4F), however, show that the network is dynamically balanced [5,6] rather than structurally balanced (Fig. 4D, red dot). Each of the critical modes exhibits a slowly decaying autocorrelation function and evolves in a specific direction within the high-dimensional space of all neurons that is different from the population activity (Fig. 4F). The experimental identification of these hidden modes presents a major challenge for the future due to the strong subsampling of simultaneous single-neuron activities.
We expect implications of the operation in the dynamically balanced, critical regime (λ max 1) for learning and information processing. Neurons belonging to a critical mode show pairwise covariances that strongly exceed  the average (Fig. 4E). Such covariances interact with spike-timing dependent synaptic plasticity [29], leading to an increased interaction between neuronal and synaptic dynamics. Weak external inputs to the network are, moreover, sufficient to shift a macroscopic number of eigenvalues across the edge of stability [16] and are therefore potent to modulate the integrative properties of the network. Critical modes have a multitude of characteristic shapes and life times (Fig. 4F inset, [30]); they here arise despite the stereotypical and fast dynamics of individual neurons as a result of the disordered network structure. This diversity is in stark contrast to the small number of different modes at the point of close structural balance (here: two modes, Fig. 4C inset). The enriched repertoire enables the parallel integration and maintenance of information over prolonged time scales. Such networks provide a wealth of transformations on the input and therefore may serve as an exhaustive reservoir for computation [31,32].
Finally, the two types of criticality are not mutually exclusive as they are governed by different mechanisms. They can coexist in different brain regions or even in the same local network; networks may hence dynamically move into either regime, adapting brain function to momentary demands.

I. MEAN-FIELD THEORY FOR META-STATISTICS BEYOND SELF-AVERAGING
Mean-field theory most-commonly employs the thermodynamic limit (N → ∞), reducing the collective dynamics of the N interacting units to N pairwise independent units, each subject to a self-consistently determined auxiliary field [4,17,33,34]. Covariances of individual neurons are self-averaging in this limit, so they are identical for all units and independent of the realization of the randomness in the network connectivity. In particular, cross-covariances vanish. In the absence of external stimuli and in the weakly correlated regime, covariances can be understood in linear response theory [35]. For such a linearized network model, interactions between neurons can be included as a finite-size correction within this self-averaging framework [5,6,20,21,36]. While this procedure yields covariances averaged over many pairs of units, the experimentally observed neuron to neuron variability [8] is lost. Here, we develop a theory beyond self-averaging covariances that captures their statistics. It exploits the large size of biologically realistic local networks, but includes finite-size fluctuations of auxiliary fields which derive from the quenched disorder of the couplings to explain the variability of covariances. To perform this qualitative step, we need to combine methods from different fields: We use a functional formulation of Gaussian processes that originates from statistical field theory [22,23,37] and combine it with methods typically used for disordered systems in the large N limit [25], such as spin glasses [24,38]. We outline these steps in detail below.

A. Moment generating functional for the network dynamics
In the following, we consider time-lag integrated covariances with the moment-generating functional Here is a purely imaginary response field, Dx, Dx are suitably defined path integral measures, andx Tx = i x i (t)x i (t) dt is a scalar product [22,23,37]. The generating functional can easily be interpreted in Fourier domain due to the linearity of Eq. (6) and the invariance of scalar products under unitary transforms with Fourier transformed variables denoted by capital letters. The scalar product in frequency domain readsX T X = i X i (−ω)X i (ω) dω. The generating functional factorizes into generating functions for each frequency ω. As shown in Eq. (5), time-lag integrated covariances only require the knowledge of X(0). In the following, we will therefore only discuss zero frequency. After integration over all non-zero frequencies one obtains the generating function for zero frequency with the single-frequency (ω = 0) scalar product defined asX T X = iX i X i , and integration measures DX = B. Self-averaging meta-statistics Equation (11) relates covariances between individual pairs of neurons to the connectivity matrix W . These, however, change between realizations of the random connectivity. In contrast, the meta-statistics, by which we denote the moments of the distribution of covariances, can be assumed constant across realizations (self-averaging, [38]). We therefore seek for an expression relating the moments of the distribution of covariances to the statistics of the connectivity W .
We denote with¯the empirical average, with x the expectation over realizations of the processes x, and with the ensemble average over the disordered connectivity W . Exchanging the order of differentiation and averaging allows expressing second moments of covariances as derivatives of a single disorder-averaged generating function Z(J) : We first assume the empirical average to be self-averaging We then use its definition as ii , and that of the second cumulant c ij = X i X j x of the zero frequency components X = X(ω = 0) = ∞ −∞ x(t) dt of Ornstein-Uhlenbeck processes x to obtain As the action Eq. (9) for a single realization of W is quadratic, Wick's theorem applies, For the special case i = j, we may identify the squared second moment (13) with a fourth moment. The latter can be expressed with the help of the generating function Wick's th.
In the last step, we exchanged the order of derivatives and the expectation value over network realizations and used the symmetry of the disorder-averaged network over units. Analogously follows for i = j and N − 1 ≈ N , For clarity, we here assume independent and identically distributed weights W ij . In the resulting cumulant expansion [38,[40][41][42] κ k is the k-th cumulant for a single connection W ij [43]. For fixed connection probability p, the number of inputs to a neuron scales with the network size N . To keep the input and its fluctuations within a certain dynamic range when increasing the network size, we require synaptic weights to scale with 1/ √ N [4,44], such that the cumulant expansion is an expansion in 1/ √ N . A truncation at the second cumulant (∝ N −1 ) maps W to a Gaussian connectivity N (µ, σ 2 /N ) so that with a homogeneous mean connection weight µ ij = µ = O(1/ √ N ). The second cumulant (σ 2 /N ) is the first nontrivial contribution to the second moment of covariances. While higher cumulants of the connectivity have an impact on higher moments of the distribution of covariances, their effect on the first two moments is suppressed by the large network size.

D. Auxiliary-field formalism
The interaction term V prevents an exact calculation of the disorder-averaged generating function. A converging perturbation series can be obtained in the auxiliary-field formulation [25], where a field Q 1 = σ 2 N X T X is introduced for the sum of a large number of statistically equivalent activity variables. Using the Hubbard-Stratonovich transformation e σ 2 2NX one obtains a free theory, which is a quadratic action in the activity (X) and response variables (X), on the background of fluctuating fields Q Z Q (J) = DX DX exp S 0 (X,X) + 1 2 Q 1X TX + Q 2 X T X + J T X .
The high dimensional integrals of the free theory Z Q (J) can be solved analytically yielding a two-dimensional interacting theory in the auxiliary fields Q 1 and Q 2 . The auxiliary field formalism translates the high-dimensional ensemble average over W to a low-dimensional average over Q; it maps the local disorder in the connections to fluctuations disorder average aux. of global fields Q interacting with a highly symmetric all-to-all connected network, illustrated in Fig. S1. Only in the special case of vanishing mean connection strength µ = 0, the system factorizes into N unconnected units, each interacting with the same set of fields Q. The all-to-all network not only captures the autocovariance of a single neuron, but also the cross-covariance with any other neuron.
and expanding Y (Q 1 , Q 2 , J) around Q * (J), yields an expansion in powers of 1/N [25,45]. We focus only on the leading-order contribution Y (Q * 1 (J), Q * 2 (J), J) which describes tree-level diagrams in the Q-theory and drop all source dependence in higher Taylor coefficients with the corresponding action Note that this Gaussian theory in X,X still contains in the latter two terms contributions from the quartic interaction term V of the original theory. Variability in Q * (J), through their source dependence, therefore gives rise to cumulants of the activity variables beyond second order. The dependence of Q * (J) on external sources J was neglected in prior work [24] since it scales to leading order as 1/N . This approximation, however, yields a Gaussian theory (see Eq. (20) for J = 0), which does not generate fourth cumulants and therefore no distribution of covariances. Taking into account the J dependence of auxiliary fields in combination with the relation between fourth cumulants and distributions of covariances Eq. (14) and Eq. (15) thus extends mean-field theory beyond self-averaging covariances.
The action (20) shows that Q * 1 (J) acts as a global contribution to the Gaussian noise whereas Q * 2 (J) directly contributes to the inverse of the covariance matrix. By integrating out the response variablesX, one can alternatively interpret both Q α * (J) as contributions to the covariance matrix of the noise. Saddle points Q * (J) are obtained self-consistently from Eq. (18) which reduces to the set of equations The above set of equations cannot be solved analytically for Q * (J). However, moments of activity result from derivatives of Z(J) evaluated at J = 0. The derivatives act not only on the source term J T X, but also on the Jdependence of the saddle-point. Therefore, moments are determined by saddle-points and their derivatives evaluated at zero source. Setting J = 0 in Eq. (21) yields the self-consistent result with where we used that , we find that second moments of activity, the average covariances, are not influenced by the source dependence of saddle points The non-zero mean connection strength µ = O(1/ √ N ) yields cross-covariances between neurons due to the finite size of the network. In addition, the fourth moments of activity are crucial for the non-vanishing variance of covariances: Formally, two of the four derivatives act on Q * (J) while the other two act on the source term J T X to yield The latter two terms in Eq. (27) correspond to the trivial Wick decomposition from the Gaussian part of the theory at Q * α = Q * α (J = 0), whereas the second derivatives of saddle points

F. Mean and variance of the covariance distribution
We obtain the mean integral covariances (see Eq. (26)) and the variance of integral covariances (see Eq. (14), Eq. (15) and Eq. (27)) The mean connectivity µ enters the coefficients γ ij = δ ij +γ (with δ ij the Kronecker symbol), χ ij = 1 N (1 + δ ij + O(1/N )) and γ = O(1/N ) only in their sub-leading corrections, and acts as a negative feedback in inhibitory or inhibitiondominated networks [6]. While this feedback suppresses mean cross-covariances (28) which consequently scale as c ij ∼ 1 N , it only yields a subleading contribution to the dispersion (29). The spread of individual cross-covariances is determined by fluctuations in connection weights and shows a scaling as δc 2 ij ∼ 1 √ N . These fluctuations, which formally originate from the variability of the auxiliary fields Q, are therefore much larger than the mean; they cause broad distributions of cross-covariances of both signs even in a homogeneous network. The expressions therefore explain the first two moments of the experimentally observed distribution of cross-covariances: mean cross-covariances scale as O(1/N ), whereas the standard deviation only scales as O(1/ √ N ). The width of the distribution is thus much larger than the mean for large networks and the distribution is centered approximately around zero.
We note that this approach is inherently different from mean-field theories for single realizations of network connectivities which keep the site dependence to infer relations between covariances and connections on the level of individual neurons [46]. In contrast, we derive a relation between the statistics of the structure and the statistics of the dynamics using a quenched average. The crucial feature of the presented theory is that heterogeneity across neurons can still be extracted from the ensemble description. We here show that the heterogeneity is closely linked to fluctuations of the auxiliary fields, which are accessible by studying their source dependence.

II. BIAS AND ERROR OF ESTIMATING THE DISPERSION IN EXPERIMENTAL DATA
Due to limited numbers of trials, the estimator for the dispersion of the covariance may be biased towards larger variances. We seek to find a correction for this bias, as illustrated in Fig. S2a. To this end we derive the dependence of experimentally estimated moments of the activities on the number of neurons and trials. The derivation here follows the standard approach of correcting the estimation, sometimes called Bessel's correction [47].
We assume a probability distribution p(n 1 , ..., n NT ), n k ∈ N N 0 , of activities n k i of neuron i in trial k ∈ {1, ..., N T }. We assume that the spike counts are sufficiently large so that their statistics is to leading order described by its first two cumulants, which is in line with our findings in Section I. We further assume that there are no correlations between different trials, so that the probability distribution factorizes over trials and the moment generating function reads where m is the vector of mean activities with mean m and variance δm 2 across neurons, c is the covariance matrix with mean autocovariance a, variance of autocovariances δa 2 , mean cross-covariance c and variance of cross-covariances δc 2 across neurons. The meta-statistics are assumed to be identical across trials. Furthermore, we assume the metastatistics to be the same for all realizations of distributions of mean activities and covariances across neurons. In the following, denotes the average over these realizations obtained from the averaged moment generating function andˆdenotes the empirical estimates of mean activities and covariances from N recorded neurons in N T trials of the experiment. Note that the average across realizations formally introduces correlations between different trials (mixed k, l terms in (31)) which allow us to calculate corrections due to the finite number of trials. Using Eq. (31), it is straight-forward to show and well known that an empirical covariance defined asĉ b ij = 1 NT NT k=1 (n k i −m i )(n k j −m j ) yields a biased estimator: withm i = 1 NT NT k=1 n k i and n k i n l j = δ ij (δ kl a + δm 2 ) + (1 − δ ij )δ kl c + m 2 . The unbiased estimator is therefore given byĉ ij = NT [47]. Using this definition, along the same lines a lengthy, but straightforward analogous calculation shows that the variance of cross-covariances (i = j) defined asδ c 2 = of the true variance of cross-covariances δc 2 . For a finite number of trials, there is a significant bias of δ c 2 caused primarily by the average variance a of spike-counts across trials (see Fig. S2a for an empirical estimate of the bias with trial-shuffled data). The aim is hence to correct for the bias of the estimator. This can be formally done using Eq. (33), which we validate for synthetic data in Fig. S2b. A fit of Eq. (33) to the variance of cross-covariances obtained from subsampled numbers of trials of the experimental data shows that the true variance of cross-covariances can be approximately obtained by subtracting the mean of the distribution of the surrogate data in Fig. S2a.
In conclusion, the finite number of trials yields a significant bias to the variance of cross-covariances, but not to the mean autocovariances, if properly defined (see Eq. (32)). Fitting data to Eq. (33) and extrapolating for N T → ∞ yields the unbiased estimate which we use in the main text. However, even neglecting this correction, the order of magnitude of the ratio ∆ between the two estimates, which determines the spectral radius (Eq. (4) in the main text) and the operational regime, is not changed by the bias.
In addition to the bias in the estimation of the variance, the estimator comes with a statistical uncertainty due to the limited number of recorded neurons. We observe n = N (N − 1) cross-covariances of which we want to estimate the true variance. The relative standard error of the variance is therefore given by 2/(n − 1) ≈ 0.009 < 1% [48]. Error propagation to the estimation of R 2 shows that the relative error of R 2 is even smaller for R 2 > 1/3, and hence in the critical regime close to R 2 ≃ 1 this error is negligible.

III. NUMERICAL SIMULATIONS
A. Figure 2 To illustrate the general mechanism relating network heterogeneity, eigenvalues, and distributions of covariances, we consider a simple network model with a sparsely and randomly connected inhibitory population of size N = 1000 and covariance of single-unit fluctuations D = 1. Inhibitory network with distance-dependent connectivity: We consider a network of N = 10000 inhibitory neurons (D = 1) randomly positioned on a 1 mm × 1 mm sheet. Each neuron receives K = 100 incoming connections of uniform strength varied in the range w = −0.1, ..., −0.01. Connections are drawn from a connection profile p(x) ∼ exp(−x 2 /(2∆ 2 )) where x is the Euclidean distance between the presynaptic and postsynaptic neuron, and ∆ = 50 µm the space constant.
Leaky integrate-and-fire network: We consider a network of N = 10000 inhibitory leaky integrate-and-fire model neurons with delta-shaped postsynaptic currents. The membrane potential of each neuron follows the differential equation where s j (t) = k δ t − t j k is the spike train of the j-th neuron and t j k denotes the k-th spike of neuron j, which occurs whenever V j exceeds the threshold θ. The membrane potential is reset V j (t j k +) ← V r to the reset potential V r after each such event and held at this level for the absolute refractory time τ r . The time lag h equals the resolution of the time-driven simulation. We choose a connection probability p = 0.1 and uniform weights varied in the range J = −1.1, ..., −0.1. Neuron and simulation parameters are shown in Table T1.
C. Figure 4 Linear instability is determined by an eigenvalue close to the critical line Re(λ) = 1. For networks with close structural balance (Fig. 4 left), this eigenvalue corresponds to the population activity and arises from a marginal net excitation of order 1/N , where N is the network size. The simplest network that shows this feature is a network of excitatory neurons. We consider a sparse, random network of N = 1000 neurons with independent and identically distributed connections W ij   Table T1. Specification of neuron and simulation parameters for network of LIF model neurons shown in Fig. 3 of the main text.
is λ max = N p(1 − p)w 2 = 0.93 1. Due to the finite size of the network, the largest real eigenvalue, which is used in Fig. 4F, slightly differs from this result and is approximately the same as for the excitatory network (λ ≈ 0.917).