Scaling of causal neural avalanches in a neutral model
SScaling of causal neural avalanches in a neutral model of neurons
Sakib Matin, ∗ Thomas Tenzin, and W. Klein
1, 21
Department of Physics, Boston University,Boston, Massachusetts 02215, USA Center for Computational Science, Boston University,Boston, Massachusetts 02215, USA (Dated: June 16, 2020)
Abstract
Neural avalanches are collective firings of neurons that exhibit emergent scale-free behavior.Understanding the nature and distribution of these avalanches is an important element in under-standing how the brain functions. We study a model of the brain for which the dynamics aregoverned by neutral theory. The neural avalanches are defined using causal connections betweenthe firing neurons. We analyze the scaling of causal neural avalanches as the critical point is ap-proached from the absorbing phase. By using cluster analysis tools from percolation theory, wecharacterize the critical properties of the neural avalanches. We identify the tuning parametersconsistent with experiments. The scaling hypothesis provides a unified explanation of the powerlaws which characterize the critical point. The critical exponents characterizing the avalanchedistributions and divergence of the response functions are consistent with the predictions of thescaling hypothesis. We use an universal scaling function for the avalanche profile to find that thefiring rate for avalanches of different sizes shows data collapse after appropriate rescaling. Wealso find data collapse for the avalanche distribution functions, which is a stronger evidence ofcriticality than just the existence of power laws. Critical slowing-down and power law relaxationof avalanches is observed as the system is tuned to its critical point. We discuss how our resultsmotivate future empirical studies of criticality in the brain. ∗ [email protected] a r X i v : . [ phy s i c s . b i o - ph ] J un . INTRODUCTION Systems with many interacting units can exhibit phenomena at macroscopic scales whichcannot be elucidated from their microscopic behavior [1, 2]. For a system at its critical point,emergent phenomena occur at all length scales and can be understood using concepts such asscaling and universality [1, 3, 4]. There is a growing interest in the question whether certainbiological systems operate near a critical point [5–9]. One question that has received muchattention is the question of whether the brain operates near a critical point [5, 10–18] and iscommonly referred to as the criticality hypothesis [5]. Experiments have shown that neuralavalanches in vivo and in vitro can exhibit scale-free behavior similar to thermal systemsnear the critical point [19–24]. The interest in the criticality hypothesis has been amplifiedby arguments that criticality in the brain may benefit memory storage and informationprocessing [5, 25–34].Although scale-free avalanches have been observed in real neural systems [20], the exis-tence of a tuning parameter is disputed [35]. An alternate explanation is that neural systemsmay exhibit self-organized criticality [36–38] where systems generically exhibits power lawdistributions without the need for any tuning [39, 40]. However, self-organized criticalityhas strict requirements such as time scale separation and local conservation laws which maynot apply to neural systems [41, 42].The temporal-proximity binning method of defining neural avalanches [20–24, 43, 44] canshow discrepancies from the true behavior of avalanches especially when multiple avalanchespropagate through the system [41, 45]. Reference [41] partly addressed this issue by intro-ducing a model where the dynamics are determined by neutral theory and the avalanchesare defined by tracking causal-connections between firing neurons. We refer to the modelintroduced in Ref. [41] as the Neutral Model of Neurons (NMN). The avalanche distribu-tions in the active phase of the NMN was studied in Ref. [41]. In this paper we study thescaling of causal neural avalanches as the the critical point is approached from the absorbingphase. Our analysis reveals there are two relevant scaling fields or tuning parameters in theabsorbing phase thereby ruling our self-organized criticality in the NMN. We show that thecausal neural avalanches in the absorbing phase are consistent with the scaling hypothesis.Additionally, we discuss how our analysis can motivate future experiments.The remainder of the paper is structured as follows. In Sec. II, we outline the neutral2odel of neurons and discuss the connections to experimental studies of criticality in neuralsystems. We provide a brief pedagogical introduction to cluster scaling methods in Sec. III.In Sec. IV, we measure the critical exponents τ and σ which characterize the scale-freecausal avalanche distributions at the critical point. We approach the critical point fromthe absorbing phase. We also find data collapse of the distributions of avalanche size andduration near the critical point. In Sec. V we study how the response function diverges as apower law as the critical point is approached. We show that our measured critical exponentsare consistent with the scaling hypothesis in Sec. VI. We find that the relevant scaling fields inthe NMN are consistent with the different tuning parameters in experiments [20, 21, 46, 47].In Sec.VII, we analyze the scaling relation between the avalanche size and duration. InSec. VIII, we use scaling arguments to derive the universal avalanche profile and show datacollapse for the firing rates for avalanches of different sizes. Critical slowing down is analyzedin Sec. IX. In Sec. X, we analyze the universal relaxation dynamics in the NMN near thecritical point. Lastly, in Sec. XI we discuss our results, which indicate that the NMN isconsistent with the predictions of the scaling hypothesis. Additionally, we discuss how ourresults may inform future empirical studies of neural systems.3 I. MODEL
The brain is a complex system consisting of many interacting neurons. In the restingstate, neurons have intrinsic voltages which fluctuate around some residual value. A neuron,triggered by some stimuli (endogenous or external), sends its action potential or spikesto its connected neighbors. A recipient neuron may also fire and send its spike to theconnected neighbors thereby resulting in an avalanche. The firing neurons in this avalancheare causally connected. Reference [41] defined neural avalanches in the neutral model ofneurons using causal-connections and showed that this definition of neural avalanches doesnot suffer the ambiguities commonly found in the temporal-proximity binning method [45].Multiple neural avalanches can propagate concurrently in the NMN because of the causallyconnected definition. The avalanches are neutral [41, 48, 49] or symmetric because the ratesthat describe the dynamics of the avalanches are the same for all labels, which distinguishthe different avalanches.The neutral model of neurons consists of N neurons which are fully connected. Everyneuron interacts with every other neuron. A neuron is either inactive, I , or active A k ,where the index k denotes the avalanche label. In Fig. 1, different colors correspond todifferent causal neural avalanches. The stochastic dynamics of the avalanches are describedby rate equations. A new avalanche with a new label is triggered at the driving rate (cid:15) . Anavalanche increases in size at the propagation rate λ as inactive neurons are triggered byactive neurons. Active neurons become inactive at the decay rate µ . The rate equationsdescribing the NMN [41] are I (cid:15) −→ A max[ k ]+1 (1) I + A k λ −→ A k + A k (2) A k µ −→ I (3)An avalanche ends when all neurons with a given label k become inactive. The size S of thecausal neural avalanche is the number of activations and the avalanche duration D is thetime between the activation of the first neuron with label k to when all neurons with index k become inactive, as shown in Fig. 1.The NMN captures many of the salient biological mechanisms relevant to neural avalanches.In neural systems, the ratio of inhibitory to excitatory neurons [20, 46, 47] and the sponta-4 IG. 1.
Neutral theory describes the dynamics of causal neural avalanches . The boxescorrespond to neurons and each column to the right is the system at a later time step. Differentcolors correspond to different causal avalanches. The rates for the dynamics are identical for allavalanches. A new avalanche is triggered at the driving rate (cid:15) and is marked with a star. An activeneuron can trigger an inactive neuron anywhere in the system at the propagation rate λ . Bothneurons share the same label as they are causally connected. An active neuron becomes inactiveat the decay rate µ . For neural avalanche labeled by black, the size S = 2 is the total number ofactivated neurons and the duration D = 2 is the time elapsed between the first activation until allthe black neurons become inactive. neous triggering rate [21, 50] affect the statistics of neural avalanches. In the NMN, we cantune the propagation rate, λ , and the decay rate, µ , to achieve a similar result to varying theratio of inhibitory and excitatory neurons in experiments. In addition, the driving rate (cid:15) inthe NMN is analogous to the spontaneous triggering of neurons [21, 51]. We can explore thecriticality of neural avalanches in the neutral model of neurons in a way that is comparableto experiments.In this paper, we only study causal neural avalanches in the absorbing phase of the NMN,where µ ≥ λ and (cid:15) ≥
0. We simulate the model using a discrete-time asynchronous randomupdate [52, 53]. All our simulations are run at N = 10 , similar to Ref. [41].5 II. THEORY
Different physical systems exhibit universal properties near the critical point [3]. We canunderstand the nature of a critical point by analyzing the statistical properties of fluctua-tions or clusters in the system [1, 54, 55]. By using real-space renormalization group tech-niques [4, 56], we can map domains in magnetic systems undergoing a thermodynamic phasetransition to clusters in percolation models undergoing a geometric phase transition. Clus-ter analysis methods from percolation theory have even been used to study non-equilibriumphase transitions in integrate-and-fire systems [57]. Reference [58] provides a summary ofscaling in percolation theory. We will use similar cluster analysis methods to study thescaling of causal avalanches in the absorbing phase of the neutral model of neurons.We can obtain thermodynamic quantities from the the avalanche number density, n S ( S )which characterizes the probability of an avalanche of size S . Near the critical point, thenumber density for S (cid:29) n S ∼ S − τ G (cid:18) SS c (cid:19) . (4)The τ exponent characterizes the power law distribution at the critical point. The char-acteristic size S c diverges as the system approaches the critical point, and the avalanchenumber density becomes a power law, n S ∼ S − τ . The Fisher ansatz [59], assumes that thescaling function G ( x ) is given by G ( x ) = exp( − x ). Equation. 4 and the Fisher Ansatz hasbeen applied to many systems [55, 57] and we will show that it is consistent with the neutralmodel of neurons.A system is said to be scale-free when the avalanche distribution follows a true power law, n S ∼ s − τ . In the scale-free case, we can plot n S = N s − τ on a log-log plot to find a straightline with slope − τ . Deviations from the power law occur when the system is not at thecritical point. The exponential cutoff corresponds to G ( x ) and appears as a “knee” in thelog-log plot [57]. Distributions with a cut-off at some characteristic scale are not scale-free.The finite lattice size in simulations also sets a cut-off.The scaling hypothesis was originally introduced to explain the universal behavior acrossdisparate thermodynamic systems [1]. According to the scaling hypothesis, the asymptoticbehavior of various thermodynamic functions follow power laws characterized by critical ex-ponents. The power laws are interconnected because they are caused by the same underlying6echanism. The predictions of the scaling hypothesis have been verified in both experimentsand numerical models [1]. Later, renormalization group methods have been used to justifythe scaling hypothesis [3].We review the scaling of thermodynamic quantities and introduce critical exponents fromstatistical mechanics. Response functions quantify the change in macroscopic behaviorcaused by changes in intensive parameters [1]. Our definition of the response function χ is the same as in percolation theory [54, 58, 60], namely χ = (cid:80) S S n S (cid:80) S S n S . (5)The response function χ is analogous to the magnetic susceptibility in the Ising model [1].The response function diverges [1] as χ ∼ h − γ , (6)where we have denoted the tuning parameter by h . Divergent response functions are hall-marks of a system near its critical point [1]. The characteristic avalanche size S c scales [1, 54]as S c ∼ h − /σ . (7)We use scaling arguments to relate the asymptotic behavior of thermodynamic quantitiesnear a critical point. The scaling of the characteristic avalanche size S c as a function thecorrelation length serves as a good pedagogical example. The correlation length ξ scales as ξ ∼ h − ν = ⇒ ξ − /ν ∼ h (8)The characteristic avalanche size scales as S c ∼ h − /σ ∼ (cid:2) ξ − /ν (cid:3) − /σ ∼ ξ /σν (9)We will show that similar scaling arguments hold for the causal avalanches in the neutralmodel of neurons. 7 V. SCALE-FREE AVALANCHE DISTRIBUTION
In this section, we study the causal avalanche size and duration distributions in theabsorbing phase of the NMN where µ ≥ λ and (cid:15) ≥ in vitro and in vivo show that thespontaneous triggering rate of neurons may be interpreted as a tuning parameter. Thisparameter corresponds to (cid:15) in the NMN. The experiments reported in Ref. [20, 46, 47] usepharmacological means to alter the excitation-inhibition ratio to alter the proximity to thecritical point. We can achieve similar results by varying the rates λ and µ . Our analysisshows that the relevant scaling fields for the NMN are the driving rate (cid:15) and propensity∆ = µ − λ , which are consistent the different experimental results [20, 21, 46, 47]. We findthat the critical point for the NMN is ∆ = µ − λ = 0 and (cid:15) = 0. We find that the distri-butions of avalanche sizes and durations follows power laws for ∆ = 0 and (cid:15) = 0. Largeravalanches are suppressed for (cid:15) > > n S and duration n D satisfypower law at the critical point with the exponents τ and τ D . As we tune the system awayfrom the critical point by increasing the driving (cid:15) or the propensity ∆, we find exponentialsuppression of the large avalanches characterized by the exponents σ and σ D for the size andduration respectively. The distribution functions for the avalanche size and duration are n S ∼ S − τ exp (cid:20) − SS c (cid:21) (10) n D ∼ D − τ D exp (cid:20) − DD c (cid:21) (11)where, S c is the characteristic avalanche size and D c is the characteristic duration size [54,58, 59]. For critical propensity, ∆ = 0, the characteristic avalanche size scales as S c ∼ (cid:15) − /σ and the characteristic duration scales as D c ∼ (cid:15) − /σ D . When (cid:15) = 0 the scaling is S c ∼ ∆ − /σ and D c ∼ ∆ − /σ D .Our measured value of τ and τ D in Fig. 2 are consistent with the theoretical mean fieldvalues, τ MF = 1 . τ MF D = 2 [5] and experimentally reported values [21]. In Fig. 3 we fitthe exponents, σ and σ D which characterize the exponential suppression of large avalanches.The critical exponents σ and σ D for neural systems have not been reported in experiments.8 remarkable consequence of the scaling hypothesis is the existence of universal scalingfunctions which are usually obtained via data collapse [1, 3, 61], where results for differentvalues of the control parameters collapse on to a single curve after appropriate rescaling.In Fig. 2, the insets show data collapse for the distributions of the causal avalanches. InFig. 2A, the causal avalanche size distribution for different values of (cid:15) is plotted. The scalingform is n S ∼ S − τ exp (cid:18) − SS c (cid:19) ∼ S − τ exp (cid:0) − S(cid:15) /σ (cid:1) (12)We multiply both sides by (cid:15) − τ/σ , n S (cid:15) − τ/σ ∼ (cid:15) − τ/σ S − τ exp (cid:0) − S(cid:15) /σ (cid:1) (13) n S (cid:15) − τ/σ ∼ (cid:0) S(cid:15) /σ (cid:1) − τ exp (cid:0) − S(cid:15) /σ (cid:1) (14)In the inset of Fig. 2A, we plot n S (cid:15) − τ/σ as a function of the rescaled avalanche size S(cid:15) /σ tofind data-collapse. The exact same scaling arguments be used to derive the rescaled scaledvariables in Fig. 2B, 2C and 2D, which are as follows n D (cid:15) − τ D /σ D ∼ (cid:0) D(cid:15) /σ D (cid:1) − τ D exp (cid:0) − D(cid:15) /σ D (cid:1) (15) n S ∆ − τ/σ ∼ (cid:0) S ∆ /σ (cid:1) − τ exp (cid:0) − S ∆ /σ (cid:1) (16) n D ∆ − τ D /σ D ∼ (cid:0) D ∆ /σ D (cid:1) − τ D exp (cid:0) − D ∆ /σ D (cid:1) (17)The data collapse of the causal avalanche distributions in the insets of Fig. 2 is compellingevidence that (cid:15) and ∆ are the relevant scaling fields or tuning parameters in the absorbingphase neutral model of neurons. Data collapse of the avalanche distributions can be used asa test of criticality in future experiments with neural systems.9 avalanche size, S p r o b . o f s i z e S , n S AAA = 0= 0.01= 0.001= 0.0001 avalanche duration, D p r o b . o f d u r a t i o n D , n D BBB = 0= 0.01= 0.001= 0.0001 avalanche size, S p r o b . o f s i z e S , n S CCC = 0= 0.01= 0.001= 0.0001 avalanche duration, D p r o b . o f d u r a t i o n D , n D DDD = 0= 0.01= 0.001= 0.0001 FIG. 2.
The distribution of avalanche sizes and durations follows a power law at thecritical point, ∆ = 0 and (cid:15) = 0 . There is exponential suppression of large avalanches for ∆ > (cid:15) >
0. The various exponents are the same for both tuning parameters. The inset showsdata-collapse for the neural avalanche distributions. In the top row , ∆ = 0 and (cid:15) is varied. (A)The exponents for the avalanche size distribution τ = 1 . ± .
05. (B) The avalanche durationdistribution is characterized by τ D = 1 . ± .
11. In the bottom row , ∆ is varied for (cid:15) = 0. (C)The corresponding exponent τ = 1 . ± .
05. (D) τ D = 1 . ± .
11 is the critical exponent for theavalanche duration distribution. Driving, C h a r a c t e r i s t i c S i z e , S c A = 0 fit 10 Driving, C h a r a c t e r i s t i c D u r a t i o n , D c B = 0 fit10 Propensity, C h a r a c t e r i s t i c S i z e , S c C = 0 fit 10 Propensity, C h a r a c t e r i s t i c D u r a t i o n , D c D = 0 fit FIG. 3.
The characteristic avalanche size S c and duration D c diverge as the criticalpoint is approached. For each value of ∆ or (cid:15) , we fit the avalanche size and duration distributionsto Eq. 10 and Eq. 11 to obtain S c and D c respectively. In the top row , (cid:15) is varied at fixed ∆ = 0.(A) The characteristic avalanche size S c diverges with critical exponent σ = 0 . ± .
07. (B)The characteristic duration D c diverges with critical exponent σ D = 0 . ± .
06. In the bottomrow , ∆ is varied at fixed (cid:15) = 0. (C) The measured exponent for S c is σ = 0 . ± .
08. (D) Thecharacteristic duration diverges with the exponent σ D = 1 . ± . . DIVERGENT RESPONSE FUNCTION Divergent response functions are hallmarks of critical points [1]. We consider responsefunction χ , which is from percolation theory, is given by [57, 58] χ = (cid:80) S S n S (cid:80) S S n S . (18)In Fig. 4, the exponent γ is found to be the same whether the critical point is approachedby decreasing ∆ at (cid:15) = 0 or by decreasing (cid:15) at ∆ = 0. When we vary ∆ we measured γ = 2 . ± .
02 and when we vary (cid:15) the exponent is γ = 1 . ± . propensity, =10 R e s p o n s e f un c t i o n , A = 0 fit 10 driving rate, R e s p o n s e F un c t i o n , B = 0 fit FIG. 4.
The response function χ diverges as the system approaches the critical point. The exponent γ characterizes the divergence and is the same when the critical point is approachedat constant driving (cid:15) = 0 or ∆ = 0. (A) For (cid:15) = 0, when ∆ is varied the exponent is γ = 2 . ± . (cid:15) , the exponent is γ = 1 . ± .
04 for ∆ = 0.
Our results indicate that the response function in the NMN diverges along either scalingfield in the absorbing phase. The requirement of tuning to the critical point implies thatthe NMN is inconsistent with self-organized criticality, for which scale-free phenomena areindependent of any tuning parameter. 12
I. SCALING RELATIONS
The scale-free behavior at critical points can be attributed to underlying singularitiesin theromodynamic functions [1, 59]. We use scaling theory to relate the different criticalexponents in the neutral model of neurons.According to the scaling hypothesis, the response functions can be described by general-ized homogeneous functions near the critical point [1]. We can write χ in Eq. 18 as χ = (cid:80) S S n S (cid:80) S Sn S (19)= (cid:82) ∞ S − τ G ( S/S c ) d S (cid:82) ∞ S − τ G ( S/S c ) d S , (20)We make change of variables u = S/S c to find χ = S − τc (cid:82) ∞ /S c u − τ G ( u ) d uS − τc (cid:82) ∞ /S c u − τ G ( u ) d u , (21)and set G ( u ) = exp( − u ) [59] to express χ as χ = S c Γ(3 − τ )Γ(2 − τ ) , (22)where Γ is the gamma function.Near the critical point, we find that the χ scales as χ ∼ S c . The scaling relation betweenthe critical exponents are (cid:15) − /σ ∼ (cid:15) − γ = ⇒ γ = 1 /σ. (23)Our measured exponents in Figs. 3 and 4 is consistent with Eq. (23). The results are the samewhen we vary ∆ at (cid:15) = 0. Our derivation is similar to scaling arguments used in percolationtheory [58], except that the denominator of the response function also contributes to thedivergence. 13 II. SIZE-DURATION SCALING OF AVALANCHES
We can relate the size of the causal avalanches to their duration using scaling argumentssimilar to Ref. [61]. The critical point is approached by varying (cid:15) for ∆ = 0. The averageavalanche size follows the same scaling as the characteristic avalanche size, (cid:104) S ( D ) (cid:105) ∼ (cid:15) − /σ .The correlation length scales as ξ ∼ (cid:15) − ν . The dynamic critical exponent z relates theduration D to ξ [62, 63]. The scaling of D is D ∼ ξ z ∼ (cid:15) − νz (24) D − /νz ∼ (cid:15) (25)Now we can relate the average size (cid:104) S ( D ) (cid:105) to the duration as (cid:104) S ( D ) (cid:105) ∼ (cid:15) − /σ ∼ (cid:2) D − /νz (cid:3) − /σ (26) ∼ D / ( σνz ) (27)The same scaling arguments can be used when ∆ is varied for (cid:15) = 0. Duration, D A v e r a g e s i z e , < S > fit FIG. 5.
The scaling of the average avalanche size as a function of the duration atthe critical point is consistent with the scaling laws.
At the critical point the scaling is (cid:104) S ( D ) (cid:105) ∼ D /σνz . The numerical estimate of σνz = 1 . ± .
03 is consistent with Eq. 28.
The critical exponents for the distribution of avalanche size and durations can be related14o the exponent σνz , by the identity [61] τ D − τ − σνz . (28)We measure the exponents on the left and right side of Eq. 28 independently. The measuredvalue of 1 / ( σνz ) = 1 . ± .
03 is consistent with ( τ D − / ( τ −
1) in Fig. 2. The scaling relationbetween the size and duration of neural avalanches have been verified experimentally [22, 23].The predictions of the scaling hypothesis provide stricter criteria for criticality than just theexistence of power law distributions. 15
III. UNIVERSAL AVALANCHE PROFILE
The avalanche profile, which describes the firing rate as a function of time, can be de-scribed by a universal scaling function near the critical point . The firing rate correspondsto the number of activations per unit time. From the scaling hypothesis, we assume the av-erage firing rate is described by a generalized homogeneous function, which can be writtenas f R ( t, D ) = D b f R ( t/D ) [1, 61, 64, 65]. We compute exponent the exponent b , (cid:104) S ( D ) (cid:105) = (cid:90) f R ( t, D ) dt = (cid:90) D b f R ( t/D ) dt ∼ D b +1 . (29)By using the scaling identity S ( D ) (cid:105) ∼ D / ( σνz ) , we find b = 1 / ( σνz ) −
1. Figure 6 shows thedata collapse for avalanches of different durations, when we scale the firing rate by D − / ( σνz ) and plot it as a function of the scaled time t/D . Our derivation follows Ref. [61].16 .0 0.2 0.4 0.6 0.8 1.0 Scaled Time, tD S c a l e d F i r i n g R a t e , f R D z D=1000D=2000D=3000
FIG. 6.
Neural avalanches have an universal avalanche profile at the critical point.
The firing rate scaled by D − / ( σνz ) as a function of the scaled time t/D shows data collapse foravalanches of different durations. Data collapse is an impressive example of universality in neural avalanches. Universalscaling functions can be used as strict criteria for criticality because the data collapse isobserved only sufficiently close to the critical point. Under certain circumstances, datacollapse for the avalanche profile has been reported for in vitro experiments [22, 23, 66],where the avalanches are defined using the temporal proximity binning method. Our resultshows that causal avalanches in the absorbing phase of the NMN follow similar scalingbehavior. 17
X. CRITICAL SLOWING DOWN
A well-known consequence of criticality is a divergent time scale [67], which is known ascritical slowing down. Here, we analyze how the time to reach a stationary state divergesas the system is tuned to the critical point.We study the equilibration time for NMN after initializing with a single active neuron.We analyze the dynamics of U ( t ), which is the number of unique causal avalanches at time t . U ( t ) is a population level quantity because it is computed using information from thewhole system at a particular instance in time. We determine that the system has reacheda stationary state when U ( t ) reaches a constant rolling time average. In Fig. 7 we plot thetime for U ( t ) to reach a steady state T E as a function of the driving rate (cid:15) at ∆ = 0. We canuse scaling arguments to relate the dynamic exponent z to the other critical exponents; thecorrelation length and time scale as ξ ∼ (cid:15) − ν and t ∼ ξ z respectively. In Fig. 7 the systemapproaches the critical point, T E diverges as T E ∼ ξ z ∼ (cid:15) − νz . (30)Our measured dynamic critical exponent νz = 0 . ± .
04 is consistent with Eq. 28 and1 / ( σνz ), where σ is determined in Fig. 5. The results remain the same when we repeat theanalysis using the total activity, ρ , instead of U . Hence, the NMN exhibits a divergent timescale characteristic of critical systems and is consistent with the scaling hypothesis.18 IG. 7.
The equilibration time T E diverges as the system approaches the critical point (cid:15) → with ∆ = 0 . T E is the the time for U ( t ) to reach a steady-state value when the NMN isinitiated with a single active neuron. We find T E ∼ (cid:15) − νz . The measured value, νz = 0 . ± . . RELAXATION DYNAMICS Power law temporal relaxation is a hallmark of critical systems [67]. We initialize theNMN with every neuron active and belonging to a unique causal avalanche and analyze therelaxation to either a fluctuating state or to an absorbing (inactive) state. In Fig. 8A, wevary (cid:15) for ∆ = 0 and analyze how the system decays to a fluctuating state. In Fig. 8B thesystem decays to an absorbing state as we have set the driving rate (cid:15) = 0. We find that thenumber of unique avalanches, U ( t ), decays as a power law for the critical value ∆ = 0, andexponentially for ∆ >
0. The critical exponent α characterizes the power law relaxation.Our measured value is α = 0 . ± .
04 in Fig. 8 and is consistent with the mean-field value α MF = 1 [52]. Inset plots in Fig. 8 show data collapse for the relaxation of U ( t ) by plottingthe rescaled variables U → U t α as a function of t → ( t | (cid:15) | ) νz and t → ( t | ∆ | ) νz respectivelyin Figs. 8A and B. This technique has been used in the study of directed percolation, whichexhibits a nonequilibrium phase transition [52, 53]. Δ = 0 ϵ = 0 FIG. 8.
The number of unique causal neural avalanches U ( t ) decays as a power law, U ( t ) ∼ t − α , at the critical point. The measured critical exponent α = 0 . ± .
04 matchesmean-field value α MF = 1. (A) For (cid:15) > U ( t ) reaches a fluctuatingstate. (B) For sub-critical propensity ∆ > (cid:15) = 0, U ( t ) decays exponentially to the absorbingstate. Inset plots show data collapse for the rescaled variables. (cid:15) the system evolves to a stationary state. The time-averaged valueof the number of unique clusters scales as (cid:104) U (cid:105) ∼ (cid:15) λ where λ = 0 . ± .
02. The divergencein the time scale to reach the stationary state was analyzed in Sec IX.By studying the dynamics of population level quantities in the NMN, we have foundcharacteristic power laws at the critical point. Furthermore, the data collapse for the decayof U ( t ) highlights the universal dynamics in the NMN near the critical point.21 I. DISCUSSION
Scale-free neural avalanches in vivo and in vitro are a remarkable emergent phenomenawhich have intrigued physicists [5] and neuroscientists [35, 68, 69]. The theory of criticalphenomena is a promising explanation of the scale-free behavior [22, 23, 70], and is furthermotivated by the arguments that criticality in the brain may have functional advantages [69].We have analyzed the scaling of causal neural avalanches in the absorbing phase of theNMN. As we allow the different parameters to approach their critical values, the responsefunction exhibit a power law divergence, similar to thermal systems [1]. Additionally, thecausal neural avalanches show scale-free distributions for ∆ = 0 and (cid:15) = 0. Large causalneural avalanches are exponentially suppressed when (cid:15) > >
0. This behavior is incontrast to systems which exhibit self-organized criticality, where power law distributionsare a generic feature. The values of the critical exponents τ and τ D in Ref. [41] are consistentwith our measured values. We used cluster analysis techniques to study the causal avalanchesand measure critical exponents, σ , σ D , γ and νz . Our measured exponents obey scaling lawsin Eq. 23. We use scaling arguments to relate the average avalanche size to the durationand numerically verify the result.We construct a strict criteria for criticality in the NMN using the scaling hypothesis. Astriking prediction of the scaling hypothesis is the existence of universal scaling functions. Inthe NMN, the avalanche profile shows data collapse after appropriate rescaling and has beenobserved in experiments [22]. We also find data collapse for the avalanche size and durationdistribution. We showed that the dynamics of population level observables in the NMN arealso consistent with criticality. When we initialize the system with a single active neuron,the steady state is reached over some characteristic time. As the critical point is approached,we find there is a divergent time scale. The dynamic critical exponent is consistent withscaling arguments. We analyze how the system relaxes from a maximally diverse state. Wefind deviations from power law decay depend on the distance from the critical point. Byusing scaling arguments, we find data collapse for the relaxation dynamics.Our analysis of the NMN shows that (cid:15) and ∆ are the relevant scaling fields in the absorbingphase. The results of Ref. [41] suggest that only (cid:15) is the relevant scaling field in the activephase of the NMN, where ∆ <
0. An important question is about the correspondencebetween the tuning parameters in experiments and those in the NMN. The experiments in22ef. [21] imply that the control parameter may be the spontaneous triggering rate, whichcorresponds to the driving rate (cid:15) in the NMN. In separate experiments [20, 46, 47], theexcitation-inhibition ratio was varied to tune the system toward criticality, which we achievedby varying the propensity, ∆ in the NMN. In experiments, the neural avalanches are definedusing temporal proximity binning method which can be different from the causal avalanchedefinition used by us and Ref. [41]. An important question is whether the scaling behaviorin the absorbing or active phase better corresponds to the experiments [20, 21, 46, 47]. Ourresults emphasize the need to incorporate causal connections in future experiments studyingneural avalanches as pointed out by Ref. [41].Our results may motivate future experimental studies of neural avalanches. The diver-gence of the response function can be used to identify the critical point in real neural systems.Additionally, we find data collapse of the avalanche profile and distributions at the criticalpoint. Experiments in Ref. [22] reported similar data collapse for avalanche profile in certainsamples.Numerous studies [5, 25, 26, 33] have discussed the possible functional benefits of criti-cality in the brain. Our results raise the important question of whether these benefits alsoapply to biological learning mechanisms such as spike-timing-dependent plasticity, whichuses causal information about firing neurons [41].Our work sheds new light on the scaling of causal neural avalanches. Our results canexplain past experimental data [20, 21], motivate questions for future studies and provide apromising path to a unified theory of neural avalanches.
ACKNOWLEDGMENTS
We would like to thank Ashish B. George and Harvey Gould for insightful comments andcritical reading of the manuscript. T.T. would like to acknowledge financial support fromBoston Universitys Undergraduate Research Opportunities Program.23
1] H Eugene Stanley.
Phase transitions and critical phenomena . Clarendon Press, Oxford, 1971.[2] Kenneth G Wilson. Problems in physics with many scales of length. 1979.[3] H Eugene Stanley. Scaling, universality, and renormalization: Three pillars of modern criticalphenomena.
Reviews of modern physics , 71(2):S358, 1999.[4] W. Klein, Harvey Gould, Natali Gulbahce, J. B. Rundle, and K. Tiampo. Structure of fluctu-ations near mean-field critical points and spinodals and its implication for physical processes.
Phys. Rev. E , 75:031114, Mar 2007.[5] Miguel A. Mu˜noz. Colloquium: Criticality and dynamical scaling in living systems.
Rev. Mod.Phys. , 90:031001, Jul 2018.[6] Thierry Mora and William Bialek. Are biological systems poised at criticality?
Journal ofStatistical Physics , 144(2):268–302, 2011.[7] Anthony A Hyman, Christoph A Weber, and Frank J¨ulicher. Liquid-liquid phase separationin biology.
Annual review of cell and developmental biology , 30:39–58, 2014.[8] Tam´as Vicsek, Andr´as Czir´ok, Eshel Ben-Jacob, Inon Cohen, and Ofer Shochet. Novel type ofphase transition in a system of self-driven particles.
Physical review letters , 75(6):1226, 1995.[9] Chung-Chuan Lo, Thomas Chou, Thomas Penzel, Thomas E Scammell, Robert E Strecker,H Eugene Stanley, and Plamen Ch Ivanov. Common scale-invariant patterns of sleep–waketransitions across mammalian species.
Proceedings of the National Academy of Sciences ,101(50):17545–17548, 2004.[10] Dante R Chialvo. Emergent complex neural dynamics.
Nature physics , 6(10):744–750, 2010.[11] L Michiels Van Kessenich, Lucilla de Arcangelis, and HJ Herrmann. Synaptic plasticity andneuronal refractory time cause scaling behaviour of neuronal avalanches.
Scientific reports ,6:32071, 2016.[12] Gerald Hahn, Adrian Ponce-Alvarez, Cyril Monier, Giacomo Benvenuti, Arvind Kumar,Fr´ed´eric Chavane, Gustavo Deco, and Yves Fr´egnac. Spontaneous cortical activity is tran-siently poised close to criticality.
PLoS computational biology , 13(5):e1005543, 2017.[13] Leonardo Dalla Porta and Mauro Copelli. Modeling neuronal avalanches and long-range tem-poral correlations at the emergence of collective oscillations: Continuously varying exponentsmimic m/eeg results.
PLOS Computational Biology , 15(4):e1006924, 2019.
14] Thomas Petermann, Tara C Thiagarajan, Mikhail A Lebedev, Miguel AL Nicolelis, Dante RChialvo, and Dietmar Plenz. Spontaneous cortical activity in awake monkeys composed ofneuronal avalanches.
Proceedings of the National Academy of Sciences , 106(37):15921–15926,2009.[15] Viola Priesemann, Michael Wibral, Mario Valderrama, Robert Pr¨opper, Michel Le Van Quyen,Theo Geisel, Jochen Triesch, Danko Nikoli´c, and Matthias HJ Munk. Spike avalanches in vivosuggest a driven, slightly subcritical brain state.
Frontiers in systems neuroscience , 8:108,2014.[16] Oren Shriki, Jeff Alstott, Frederick Carver, Tom Holroyd, Richard NA Henson, Marie L Smith,Richard Coppola, Edward Bullmore, and Dietmar Plenz. Neuronal avalanches in the restingmeg of the human brain.
Journal of Neuroscience , 33(16):7079–7090, 2013.[17] Silvia Scarpetta and Antonio de Candia. Alternation of up and down states at a dynamicalphase-transition of a neural network with spatiotemporal attractors.
Frontiers in systemsneuroscience , 8:88, 2014.[18] Ariel Haimovici, Enzo Tagliazucchi, Pablo Balenzuela, and Dante R. Chialvo. Brain organi-zation into resting state networks emerges at criticality on a model of the human connectome.
Phys. Rev. Lett. , 110:178101, Apr 2013.[19] C-C Lo, LA Nunes Amaral, Shlomo Havlin, P Ch Ivanov, Thomas Penzel, J-H Peter, andH Eugene Stanley. Dynamics of sleep-wake transitions during sleep.
EPL (Europhysics Let-ters) , 57(5):625, 2002.[20] John M Beggs and Dietmar Plenz. Neuronal avalanches in neocortical circuits.
Journal ofneuroscience , 23(35):11167–11177, 2003.[21] Johannes Zierenberg, Jens Wilting, and Viola Priesemann. Homeostatic plasticity and externalinput shape neural network dynamics.
Phys. Rev. X , 8:031018, Jul 2018.[22] Nir Friedman, Shinya Ito, Braden AW Brinkman, Masanori Shimono, RE Lee DeVille, Karin ADahmen, John M Beggs, and Thomas C Butler. Universal critical dynamics in high resolutionneuronal avalanche data.
Physical review letters , 108(20):208102, 2012.[23] Antonio J. Fontenele, Nivaldo A. P. de Vasconcelos, Tha´ıs Feliciano, Leandro A. A. Aguiar,Carina Soares-Cunha, B´arbara Coimbra, Leonardo Dalla Porta, Sidarta Ribeiro, Ana Jo˜aoRodrigues, Nuno Sousa, Pedro V. Carelli, and Mauro Copelli. Criticality between corticalstates.
Phys. Rev. Lett. , 122:208101, May 2019.
24] Elizabeth A Clement, Alby Richard, Megan Thwaites, Jonathan Ailon, Steven Peters, andClayton T Dickson. Cyclic and sleep-like spontaneous alternations of brain state under ure-thane anaesthesia.
PloS one , 3(4):e2004, 2008.[25] Chris G Langton. Computation at the edge of chaos: phase transitions and emergent compu-tation.
Physica D: Nonlinear Phenomena , 42(1-3):12–37, 1990.[26] Nils Bertschinger and Thomas Natschl¨ager. Real-time computation at the edge of chaos inrecurrent neural networks.
Neural computation , 16(7):1413–1436, 2004.[27] Ruedi Stoop and Florian Gomez. Auditory power-law activation avalanches exhibit a funda-mental computational ground state.
Phys. Rev. Lett. , 117:038102, Jul 2016.[28] Christopher J Honey, Thomas Thesen, Tobias H Donner, Lauren J Silbert, Chad E Carlson,Orrin Devinsky, Werner K Doyle, Nava Rubin, David J Heeger, and Uri Hasson. Slow corticaldynamics and the accumulation of information over long timescales.
Neuron , 76(2):423–434,2012.[29] Shree Hari Gautam, Thanh T Hoang, Kylie McClanahan, Stephen K Grady, and Woodrow LShew. Maximizing sensory dynamic range by tuning the cortical state to criticality.
PLoScomputational biology , 11(12):e1004576, 2015.[30] Paolo Massobrio, Lucilla de Arcangelis, Valentina Pasquale, Henrik J Jensen, and DietmarPlenz. Criticality as a signature of healthy neural systems.
Frontiers in systems neuroscience ,9:22, 2015.[31] Clayton Haldeman and John M. Beggs. Critical branching captures activity in living neuralnetworks and maximizes the number of metastable states.
Phys. Rev. Lett. , 94:058101, Feb2005.[32] Woodrow L Shew, Wesley P Clawson, Jeff Pobst, Yahya Karimipanah, Nathaniel C Wright,and Ralf Wessel. Adaptation to sensory input tunes visual cortex to criticality.
Nature Physics ,11(8):659, 2015.[33] Osame Kinouchi and Mauro Copelli. Optimal dynamical range of excitable networks at criti-cality.
Nature physics , 2(5):348, 2006.[34] Joschka Boedecker, Oliver Obst, Joseph T Lizier, N Michael Mayer, and Minoru Asada.Information processing in echo state networks at the edge of chaos.
Theory in Biosciences ,131(3):205–213, 2012.[35] John M Beggs and Nicholas Timme. Being critical of criticality in the brain.
Frontiers in hysiology , 3:163, 2012.[36] Anirban Das and Anna Levina. Critical neuronal models with relaxed timescale separation. Phys. Rev. X , 9:021062, Jun 2019.[37] Daniel Millman, Stefan Mihalas, Alfredo Kirkwood, and Ernst Niebur. Self-organized critical-ity occurs in non-conservative neuronal networks during upstates.
Nature physics , 6(10):801,2010.[38] C. B´edard, H. Kr¨oger, and A. Destexhe. Does the 1 /f frequency scaling of brain signals reflectself-organized critical states? Phys. Rev. Lett. , 97:118102, Sep 2006.[39] Per Bak, Chao Tang, and Kurt Wiesenfeld. Self-organized criticality: An explanation of the1/f noise.
Phys. Rev. Lett. , 59:381–384, Jul 1987.[40] Mark EJ Newman. Power laws, pareto distributions and zipf’s law.
Contemporary physics ,46(5):323–351, 2005.[41] Matteo Martinello, Jorge Hidalgo, Amos Maritan, Serena di Santo, Dietmar Plenz, andMiguel A. Mu˜noz. Neutral theory and scale-free neural dynamics.
Phys. Rev. X , 7:041071,Dec 2017.[42] Christof Koch and Idan Segev.
Methods in neuronal modeling: from ions to networks . MITpress, 1998.[43] Tiago L Ribeiro, Sidarta Ribeiro, Hindiael Belchior, F´abio Caixeta, and Mauro Copelli. Un-dersampled critical branching processes on small-world and random networks fail to reproducethe statistics of spike avalanches.
PloS one , 9(4):e94992, 2014.[44] Anna Levina and Viola Priesemann. Subsampling scaling.
Nature communications , 8:15140,2017.[45] Pablo Villegas, Serena di Santo, Raffaella Burioni, and Miguel A Mu˜noz. Timeseries thresh-olding and the definition of avalanche size. arXiv preprint arXiv:1902.10465 , 2019.[46] Woodrow L Shew, Hongdian Yang, Thomas Petermann, Rajarshi Roy, and Dietmar Plenz.Neuronal avalanches imply maximum dynamic range in cortical networks at criticality.
Journalof neuroscience , 29(49):15595–15600, 2009.[47] Simon-Shlomo Poil, Richard Hardstone, Huibert D Mansvelder, and Klaus Linkenkaer-Hansen.Critical-state dynamics of avalanches and oscillations jointly emerge from balanced excita-tion/inhibition in neuronal networks.
Journal of Neuroscience , 32(29):9817–9823, 2012.[48] Motoo Kimura. The neutral theory of molecular evolution: A review of recent evidence.
The apanese Journal of Genetics , 66(4):367–386, 1991.[49] Sandro Azaele, Samir Suweis, Jacopo Grilli, Igor Volkov, Jayanth R. Banavar, and AmosMaritan. Statistical mechanics of ecological systems: Neutral theory and beyond. Rev. Mod.Phys. , 88:035003, Jul 2016.[50] Viola Priesemann, Mario Valderrama, Michael Wibral, and Michel Le Van Quyen. Neuronalavalanches differ from wakefulness to deep sleep–evidence from intracranial depth recordingsin humans.
PLoS computational biology , 9(3):e1002985, 2013.[51] Christian Meisel, Eckehard Olbrich, Oren Shriki, and Peter Achermann. Fading signaturesof critical brain dynamics during sustained wakefulness in humans.
Journal of Neuroscience ,33(44):17363–17372, 2013.[52] Haye Hinrichsen. Non-equilibrium critical phenomena and phase transitions into absorbingstates.
Advances in Physics , 49(7):815–958, 2000.[53] Malte Henkel, Haye Hinrichsen, Sven L¨ubeck, and Michel Pleimling.
Non-equilibrium phasetransitions , volume 1. Springer, 2008.[54] Dietrich Stauffer and Ammon Aharony.
Introduction to percolation theory . CRC press, 2018.[55] C. A. Serino, K. F. Tiampo, and W. Klein. New approach to gutenberg-richter scaling.
Phys.Rev. Lett. , 106:108501, Mar 2011.[56] A Coniglio and W Klein. Clusters and ising critical droplets: a renormalisation group ap-proach.
Journal of Physics A: Mathematical and General , 13(8):2775, 1980.[57] Sakib Matin, Chon-Kit Pun, Harvey Gould, and W. Klein. Effective ergodicity breaking phasetransition in a driven-dissipative system.
Phys. Rev. E , 101:022103, Feb 2020.[58] Dietrich Stauffer. Scaling theory of percolation clusters.
Physics Reports , 54(1):1–74, 1979.[59] Michael E. Fisher. The theory of condensation and the critical point.
Physics , 3:255–283,1967.[60] M Sahini and M Sahimi.
Applications of percolation theory . CRC Press, 1994.[61] James P Sethna, Karin A Dahmen, and Christopher R Myers. Crackling noise.
Nature ,410(6825):242, 2001.[62] Miguel A. Mu˜noz, Ronald , Alessandro Vespignani, and Stefano Zapperi. Avalanche andspreading exponents in systems with absorbing states.
Phys. Rev. E , 59:6175–6179, May1999.[63] Andrea Roli, Marco Villani, Alessandro Filisetti, and Roberto Serra. Dynamical criticality: verview and open questions. Journal of Systems Science and Complexity , 31(3):647–663,2018.[64] James P Gleeson and Rick Durrett. Temporal profiles of avalanches on networks.
Naturecommunications , 8(1):1227, 2017.[65] Adri´an Ponce-Alvarez, Adrien Jouary, Martin Privat, Gustavo Deco, and Germ´an Sumbre.Whole-brain neuronal activity displays crackling noise dynamics.
Neuron , 100(6):1446–1459,2018.[66] Shan Yu, Hongdian Yang, Oren Shriki, and Dietmar Plenz. Universal organization of restingbrain activity at the thermodynamic critical point.
Frontiers in systems neuroscience , 7:42,2013.[67] P. C. Hohenberg and B. I. Halperin. Theory of dynamic critical phenomena.
Rev. Mod. Phys. ,49:435–479, Jul 1977.[68] Wesley P Clawson, Nathaniel C Wright, Ralf Wessel, and Woodrow L Shew. Adaptationtowards scale-free dynamics improves cortical stimulus discrimination at the cost of reduceddetection.
PLoS computational biology , 13(5):e1005574, 2017.[69] Woodrow L Shew and Dietmar Plenz. The functional benefits of criticality in the cortex.
Theneuroscientist , 19(1):88–100, 2013.[70] Anna Levina, J. Michael Herrmann, and Theo Geisel. Phase transitions towards criticality ina neural system with adaptive interactions.
Phys. Rev. Lett. , 102:118110, Mar 2009., 102:118110, Mar 2009.