Measuring neuronal avalanches in disordered systems with absorbing states
??APS/123-QED?
Measuring neuronal avalanches in disordered systems withabsorbing states
M. Girardi-Schappo
1, 2 and M. H. R. Tragtenberg ∗ Neuroimaging of Epilepsy Laboratory,McConnell Brain Imaging Center, McGill University,Montreal Neurological Institute and Hospital,H3A 2B4, Montreal, Quebec, Canada Departamento de F´ısica, Universidade Federal de Santa Catarina,88040-900, Florian´opolis, Santa Catarina, Brazil (Dated: November 5, 2018)
Abstract
Power-law shaped avalanche size distributions are widely used to probe for critical behavior inmany different systems, particularly in neural networks. The definition of avalanche is ambiguous.Usually, theoretical avalanches are defined as the activity between a stimulus and the relaxation toan inactive absorbing state. On the other hand, experimental neuronal avalanches are defined bythe activity between consecutive silent states. We claim that the latter definition may be extendedto some theoretical models in order to characterize their power-law avalanches and critical behavior.We study a system in which the separation of driving and relaxation time scales emerges from itsstructure. We apply both definitions of avalanche to our model. Both yield power-law distributedavalanches that scale with system size in the critical point as expected. Nevertheless, we findrestricted power-law distributed avalanches outside of the critical region within the experimentalprocedure, which is not expected by the standard theoretical definition. We remark that theseresults are dependent on the model details.
PACS numbers: 05.70.Jk,45.70.Ht,87.18.Hf,87.19.lj ∗ Electronic address: [email protected] a r X i v : . [ phy s i c s . b i o - ph ] M a y . INTRODUCTION Power laws (PL) are ubiquitous in nature [1], appearing in the description of phenomenaas diverse as sand pile avalanches [2, 3], rice pile avalanches [4], forest fires [5, 6], the velocitydistribution of ants in ant farms [7], the magnitude distribution of earthquakes [8], surfacegrowth [2], superconductor vortices [2], solar flares [9], stock market prices [10], connectivityof social networks [11], evolution, ecosystems, epidemics [12], neuronal avalanches [13–15],etc. Critical phenomena are characterized by PL in two different ways. In the critical point,the correlation function decays as a PL. Close to a critical point, the order parameter,susceptibility and other thermodynamical quantities have PL behavior, governed by criticalexponents [16]. Inspired by the theory of critical phenomena [16], Bak et al. [17] proposedthat these phenomena’s events are PL distributed because such systems self-organize toa non-equilibrium critical point. This was the inception of the theory of self-organizedcriticality (SOC) [2, 3].The authors used a sand pile-like cellular automaton model to describe the dynamicsof such systems: at each time step a grain of sand is added to a random pile inside asquare lattice. When a pile reaches a predefined critical slope, it redistributes four of itsexceeding grains of sand to its four neighbors; if any of the neighbors becomes unstable, italso redistributes its grains and so on, interactively, until the system reaches a new absorbingconfiguration. The result of this branching process is an avalanche of size ¯ s (total numberof redistributed grains) and duration of ¯ T time steps. As long as there is an avalanchehappening, the external deposition of grains into the system is paused. It restarts onlyafter the avalanche has ceased. This procedure is known as the separation of time scales ofexternal driving and internal relaxation. Avalanches thus generated have sizes and durationdistributed according to PL: P (¯ s ) ∼ ¯ s − ¯ α and P ( ¯ T ) ∼ ¯ T − ¯ τ , where size and duration areexpected to satisfy ¯ s ∼ ¯ T a such that ¯ α and ¯ τ follow Sethna’s scaling relation a = (¯ τ − / ( ¯ α −
1) [18].In the critical point, even in the presence of short range interactions, an initial stimuluspropagates throughout the system because of the PL long-range correlation in space andtime. It produces a 1 /f power spectrum signature in many of these systems [2]. Theexternal separation of time scales breaks down the temporal correlations between subsequentavalanches, so it is not surprising that the sand pile model generates 1 /f power spectrum2nstead [19]. The systems which have sand pile-like dynamics are known to pertain touniversality classes of phase transitions into absorbing states [3].SOC has been proposed to describe the brain activity [20, 21] and the first PL distributedneuronal avalanches were found by Beggs and Plenz [22] using in vitro local field potentials(LFP) recordings. Many models were proposed to generate avalanches using brain-inspiredsystems [23–30]. Usually, models rely on probabilistic cellular automata, dynamical maps, orintegrate-and-fire neurons connected through dynamical or stochastic synapses. All of thesemodels, however, are similar to the sand pile model in the sense that they independently gen-erate avalanches, forcefully separating stimulation and relaxation time scales. Consequently,the long-range correlations [31, 32] and 1 /f power spectrum [31, 33–36] found experimentallyfor the healthy state of the brain are not well described by most of these models. Notice,however, that some models on these classes may exhibit 1 /f -like power spectrum close tothe critical state by the introduction of inhibitory synapses [24, 37].The experimental procedure for measuring neuronal avalanches consists of low or high-pass filtering each LFP time series, and then thresholding them in order to obtain a timeseries of spike-like activity. The data are then binned in time using the average intereventinterval (each event is a spike-like activity in any electrode processed signal). Empty timebins are used to separate subsequent avalanches. Notice that these empty time bins arenot analogous to the inactive absorbing state used to separate avalanches in each of theneuronal avalanches models listed above due to the ignored subthreshold background activityof the electrodes. Such method has been applied to many works probing for neuronalavalanches [27, 30, 38–41]. In fact, experimental evidence suggests that a better analogousto such LFP recordings would be a sand pile subjected to random background activity dueto almost absence of empty time bins [41].A reasonable extension of the experimental approach can be applied to theoretical mod-els in order to probe for power-law distributed avalanches and the critical state [42–44].However, the two approaches may yield different results [44]. In this work, we present amodel in which both approaches agree in the determination of the critical state. Moreover,the second approach allows the calculation of the power spectrum for our particular modeland seems to yield restricted power-law avalanches even outside of the critical range. Noticethat these results are dependent of model parameters.Here we study a biologically plausible mammalian visual system model [45]. It presents an3bsorbing phase transition to a percolating phase in which the activity gets spontaneouslyseparated in time due to its intrinsic dynamics. An ambiguity concerning the definitionof avalanches is naturally raised: should avalanche be defined as the activity between thestimulus and the moment the system reaches the inactive absorbing state or as the networkactivity between consecutive silent states? We compare both approaches and show thatinside the critical range they are very similar and both display PL distributed avalancheswith diverging cutoff. However, outside the critical range the second approach still yieldsrestricted PL distributed avalanches whereas the first approach yields the expected breakdown of PL behavior. In fact, we show that inside the critical range, the power spectrumof the first approach is white-noise-like, while the second approach presents approximately1 /f noise [45]. II. MODEL
The model has two components: structure and dynamics. The structure is based on theform recognition pathway of the mammalian visual system [45]. There are six square layerssequentially stacked one on top of the other (Fig. 1A and B). The top layer is the inputwhere each of its elements represents a photoreceptor of the retina. The bottom layer is theoutput of the primary visual cortex (V1) where each element represents one axonal terminalthat would connect to secondary visual area of the cortex. The input and output layershave N io = (10 L ) elements each. The four internal layers are the LGN from the thalamus,and the VI, IVC β and II/III from V1, respectively. Each of these four layers contains L neurons. The network thus has N = 4 L neurons. Each neuron is composed of 100 dendriticcompartments, 1 soma compartment and 10 axonal compartments (Fig. 1C).The network is built according to the following algorithm: for each neuron of the network,(a) a postsynaptic neuron is drawn from an adjacent layer (according to arrows in Fig. 1A)inside a limited excitatory field of 7 × L L G N . . .. . . ... V I I V C β I I / I I I
InputLGNVIIVC β II/IIIOutput L L L L (10 L ) (10 L ) ABC
Layer A
Attenuation
Layer BNeuron j Neuron i ʃ ... ʃ ...
49 50 51 ...
99 101 ...
10 20 30 40 50 60 70 80 90 100
Dendritic compartment, d m P ( d m )
1 2 3 4 5 6 7 8 9 10
Axonal compartment, a k P ( a k ) (presynaptic) (postsynaptic)Synapses FIG. 1: Elements of the V1 model. A: Architecture of the network. B: Spatial organization ofthe network of N = 4 L neurons. The columnar structure is highlighted in between dotted lines.There is a column of size N c = 4 l = 196 neurons centered on each neuron of the network. C: Thecompartmental structure of the neuron: dendritic compartments (light gray), axonal compartments(black) and soma (dark gray). The synapses are formed between presynaptic axonal compartments[chosen according to an exponential probability distribution P ( a k )] and dendritic compartments[chosen with Gaussian probability distribution P ( d m )]. The signal attenuates while moving forwardin the dendrites. the center of the dendrite with standard deviation of 20 dendritic compartments (Fig. 1C);(d) the chosen axonal compartment is then connected to the chosen dendritic compartment.Each element of the input layer is connected to a dendritic compartment a neuron of theLGN in the same way that a presynaptic axonal compartment would be. The amount ofsynapses for the adjacent layers are presented in Table I. Notice that the limited excitatoryfield of presynaptic neurons in step (a) generates the well known columnar structure in the5etwork [46], illustrated in Fig. 1. One trial of the simulation consists of: building thenetwork, stimulating the network and waiting for the activity to cease. This proceduregives rise to quenched disorder. Thus, many trials of the network disorder are necessary tocalculate observables.Initially, a 30 ×
30 square of photoreceptors in the input layer is instantly activated,corresponding to approximately 3 × E , from all its active presynaptic axonal compartments with the signal from theprevious dendritic compartment. This sum is attenuated at each time step by a constant λ .A soma compartment fires an action potential if the signal in the last dendritic compartmentis greater than a threshold v T . The axonal compartments propagate the signal coming fromthe soma without any dissipation. These rules are expressed as: d ( i ) k ( t + 1) = λ (cid:34) d ( i ) k − ( t ) + E (cid:88) j,n a ( j ) n ( t ) (cid:35) , (1) v i ( t + 1) = Θ (cid:16) d ( i )100 ( t ) − v T (cid:17) , if v i ( t ) = 0 , − V R R , if v i ( t ) = 1 , v i ( t ) + 1 , if v i ( t ) < , (2) TABLE I: Amount of connections of each layer of the model per presynaptic element. Connectionscome from layers listed as rows and go to layers listed as columns of the table. Values are derivedfrom many experimental works [45].
From \ To LGN VI IVC β II/III OutputInput
LGN - - 500 - - VI - - 1100 350 - IVC β - 600 - 700 - II/III - - - - 100 = 1 : E = 1 : E = 13 : Time steps A c t i v i t y , A ( t ) E = 1 : FIG. 2: Sum of all network spikes, A ( t ), for each time step t and three typical values of E . Inset:detail of the network activity for E = 1 . a ( i ) k ( t + 1) = Θ ( v i ( t )) δ k, + a ( i ) k − ( t ) , (3)where t is the time step, d ( i ) k ( t ) [ a ( i ) k ( t )] is the potential in the dendritic [axonal] compartment k of the neuron i at time t , v i ( t ) is the soma potential of neuron i at time t , E > n of presynaptic neuron j , λ = 0 .
996 is themembrane attenuation constant, v T = 10 mV is the firing threshold, R is the refractoryperiod, V R ≡ x ) isthe Heaviside step function and δ m,n is the Kronecker delta. Notice that d ( t ) = a ( t ) ≡ E ), dissipation (controlled by λ ) and therefractory period R determines whether the system is in an inactive, active or percolatingphase [45]. R is adjusted such that there is no self-sustained activity in the interlayer loopbetween layers VI and IVC β . λ is fitted to the experimental value obtained in Ref. [47].Under these conditions, E controls a phase transition between an inactive absorbing state(subcritical activity, E < .
11 mV) and an inactive percolating phase (supercritical activity,
E > .
19 mV). The range 1 . ≤ E ≤ .
19 mV is an extended region of critical behavior [45].
III. RESULTS
The network activity for a given simulation trial, A ( t ), is simply the sum of every somaspike at each time step t , A ( t ) = N (cid:88) i =1 δ v i ( t ) , , (4)7here δ v i ( t ) , = 1 if v i ( t ) = 1 (zero otherwise) is the Kronecker delta. The profile of A ( t ) inFig. 2 allows us to define avalanches in two different ways: (I) an avalanche is regarded as allthe activity due to a single stimulus, similarly to the procedure applied to the sand pile modeland other absorbing state systems [3, 17]. In this approach, each curve in Fig. 2 for a typical E represents a single avalanche. In order to generate other avalanches, we must performother simulation trials (i.e. rebuilding the network and re-stimulating the Input layer). (II)an avalanche is all the activity surrounded by silent periods, similarly to the experimentalprocedure used to measure neuronal avalanches. Within this second approach, each curvein Fig. 2 for a typical E contains many avalanches because each avalanche corresponds toan activity peak (which are highlighted in the inset of this figure for E = 1 . E = 1 . E = 1 .
185 mV; andsupercritical E = 13 mV) in order to study the shape of avalanche distributions, P , andcomplementary cumulative distributions, F , of both avalanche sizes, s , and duration, T .The avalanche sizes and duration distributions are expected to scale according to [3]: P ( s ) ∼ s − α G s ( s/s c ) , (5) P ( T ) ∼ T − τ G T ( T /T c ) , (6)where α and τ are scaling exponents, G s,T ( x ) are scaling functions and s c and T c are thecutoff of the distributions. In the critical point, both cutoff are expected to scale with systemsize: s c ∼ L D , (7) T c ∼ L µ , (8)where L is the characteristic system size, D ( µ ) is a characteristic dimension of the avalanchesizes (duration) [3]. In fact, the PL behavior holds only for s < s c (or T < T c ) in Eqs. (5)and (6).The complementary cumulative distribution (henceforth referred only as cumulative dis-tribution) is defined as F ( s ) = (cid:82) ∞ s P ( s (cid:48) ) ds (cid:48) . In this work we rather choose to analyze thecumulative distribution because it provides a clearer visualization of the data, since it is acontinuous function of its variables. It also displays reduced noise because its precision doesnot depend on the size of the bins of the distributions histogram and it has a better defined8 valanche size, s C u m u l a t i v e d i s t ., F ( s ) -1 T F ( T ) -1 E = 1 : E = 1 : E = 13mV Avalanche size, s D i s t r i bu t i o n , P ( s ) -2 -1 E = 1 : E = 1 : E = 13mV , = 1 : T s a = 1 : BA FIG. 3: Distributions for the absorbing state avalanches (approach I) in three E regimes: ( (cid:72) )subcritical, ( (cid:108) ) critical, and ( (cid:78) ) supercritical. A: P (¯ s ) and PL fit for the critical regime yielding¯ α = 1 .
38. A inset: Sethna’s scaling relation ¯ s ∼ ¯ T a yields a = 1 . F (¯ s ) and F ( ¯ T ) for the three considered regimes. Notice that only the critical regime has PLdistributed avalanches. cutoff [48]. We may write F ( s ) and F ( T ) using scaling functions, F ( s ) ∼ s − α +1 H s ( s/s c ) , (9) F ( T ) ∼ T − τ +1 H T ( T /T c ) , (10)where H s,T ( x ) are scaling functions and s c and T c are defined in Eqs. (7) and (8).The exponents α , τ , D and µ may be calculated through the collapse of the distribu-tions [3]: we run trials for many different system sizes, L , and then plot s α − F ( s ) versus s/L D ; the best value of α collapses the data due to different L over a single horizontal line,whilst the best value of D collapses the cutoff of the distributions due to different L over asingle vertical line. This procedure may also be applied to F ( T ).Another way to calculate these exponents is to use the fact that the PL holds only for9 < s c and write F ( s ) as [28] F ( s ) = (cid:90) s c s P ( s (cid:48) ) ds (cid:48) = c + c s − α +1 , (11)where the cutoff may be directly calculated by fitting c , c and α to the computational datausing Eq. (11): s c = ( − c /c ) / ( − α +1) . The same equation may be fitted for F ( T ) yielding τ and T c . We will use a bar over s , T , α and τ to denote measurements using approach I. Avalanche duration, T
500 1000 1500 2000 2500 3000 C u m u l a t i v e d i s t ., F ( T ) T =L
10 20 30 40 T = ! F ( T ) L = 20 L = 40 L = 80 L = 99 T ! = +1 E = 1 :
185 mV = = 1 :
41= 0 : Avalanche size, s C u m u l a t i v e d i s t ., F ( s ) -1 s=L D -1 s , ! F ( s ) L = 20 L = 40 L = 80 L = 99 s ! , +1 E = 1 :
185 mV , = 1 : D = 1 : BA FIG. 4: Cumulative distributions of absorbing state avalanches (approach I) for the critical regime( E = 1 .
185 mV) and four network sizes: ( (cid:72) ) L = 20, ( (cid:108) ) L = 40, ( (cid:4) ) L = 80 and ( (cid:78) ) L = 99. A: F (¯ s ). A inset: plot of H (¯ s/L ¯ D ) ∼ ¯ s ¯ α − F (¯ s ) (collapse of the distributions) yielding ¯ α = 1 .
38 and¯ D = 1 .
75. B: F ( ¯ T ). B inset: plot of H ( ¯ T /L ¯ µ ) ∼ ¯ T ¯ τ − F ( ¯ T ) (collapse of the distributions) yielding¯ τ = 1 .
41 and ¯ µ = 0 . . Absorbing state avalanches We define the absorbing state (approach I) avalanche size, ¯ s , as the sum of all the spikesgenerated by a single stimulus (i.e. by a single simulation trial):¯ s = ¯ T (cid:88) t =0 A ( t ) , (12)where ¯ T is the avalanche duration marked by arrows in Fig. 2. ¯ T is the time neededfor the system to relax into an absorbing inactive state. Many trials of the network areconsidered, each one yields an avalanche with size ¯ s and duration ¯ T . Fig. 3A (Fig. 3B)shows the distribution (cumulative distribution) of avalanche sizes (sizes and duration) forthree different E : upside-down green triangles correspond to the subcritical phase, red circlescorrespond to the critical phase and blue triangles correspond to the supercritical phase.The PL distributions, P (¯ s ) ∼ ¯ s − ¯ α and P ( ¯ T ) ∼ ¯ T − ¯ τ , from Eqs. (5) and (6), are obtainedonly inside the critical range 1 . ≤ E ≤ .
19 mV as expected. The inset of Fig. 3A showsthat the scaling relation ¯ s ∼ ¯ T a with a = 1 . P ( s ) and P ( T ) are PL shaped. Fig. 4A and B show thecumulative distributions of avalanche size and duration for many system sizes in the criticalrange for E = 1 .
185 mV. The collapse of the cumulative distributions yields ¯ α = 1 .
38 and¯ D = 1 .
75 for the avalanche sizes, and ¯ τ = 1 .
41 and ¯ µ = 0 .
94 for the avalanche duration[Eqs. (7) to (10)]. Sethna’s scaling relation [18], a = (¯ τ − / ( ¯ α −
1) = 0 . / . ≈ . B. Silent periods avalanches
The silent periods (approach II) avalanche size, s , is defined as the sum of all the spikesof the network between two consecutive periods of silence: s = t n +1 (cid:88) t n A ( t ) , (13)where the t n are such that A ( t n ) = A ( t n +1 ) = 0 and A ( t n < t < t n +1 ) >
0. The avalancheduration is simply the number of time steps in between two consecutive moments of silence, T = t n +1 − t n [45, 49]. 11 valanche size, s C u m u l a t i v e d i s t ., F ( s ) -3 -2 -1 L
20 40 60 80 100 s c ( ) estimated s c s c L : L = 20 L = 40 L = 80 L = 99 a + bs ! , +1 B E = 1 :
185 mV , = 1 : s A v a l a n c h e s i ze d i s t ., P ( s ) -6 -4 -2 E = 1 : E = 1 :
185 mV E = 13 mV s P ( s ) -4 -2 L = 20 L = 40 L = 80 L = 99 E = 1 :
185 mV A L = 99 FIG. 5: Silent periods avalanche size distributions (approach II). A: P ( s ) for E in three regimes:( (cid:72) ) subcritical, ( (cid:108) ) critical, and ( (cid:78) ) supercritical. Notice that the three regimes have PL shapeinside the shaded area. A inset: critical regime P ( s ) for critical E = 1 .
185 mV and different L .B: F ( s ) for critical E = 1 .
185 mV and different L used to fit Eq. (11) yielding s c ( L ) and α = 1 . L = 20 (R-Squared: 99 . s c ∼ L D yielding D = 0 . The interspersing peaks of activity emerge in this model due to the deterministic dynamicsand compartmental body of the neurons. The apparent inactive moments happen whilethe signal is being propagated through dendrites and axons. The duration of such silentperiods depends on network architecture and neurons parameters, such as the amount ofdendritic and axonal compartments and the balance between the refractory period, R , andthe excitation-dissipation rate (controlled by E and λ ). The separation of avalanches is thenvery natural because there is no need to reset the system in order to spark another waveof activity. Similarly, the ignored background activity in experimental measurements (andhence the experimental separation of time scales) contains weak electric signals generatedby the propagation of potentials through membranes and synapses.12n this work, we are interested in comparing the avalanches thus defined with thoseof approach I. Notice that the avalanche sizes and duration under approach II have beenthoroughly studied and shown to be PL distributed in the critical region of this model(see Figs. 5 and 6) [45, 49]. The cutoff of such distributions, calculated here using the fitof Eq. (11) to different values of L , also non-trivially scale with system size, as expectedfor regular absorbing state avalanches (inset of Figs. 5B and 6B) [45]. However, smallavalanches are also PL distributed for sub and supercritical values of E (see the shadedarea of Figs. 5A and 6A) [45]. The distributions also obey Sethna’s scaling relation for allregimes of E [49]. This result contrasts with the distributions obtained for absorbing stateavalanches, in which only critical state avalanches display PL shape, non-trivial finite size Avalanche duration, T C u m u l a t i v e d i s t ., F ( T ) -3 -2 -1 L
20 40 60 80 100 T c ( ) estimated T c T c exp( cL ) L = 20 L = 40 L = 80 L = 99 a + bT ! = +1 E = 1 :
185 mV = = 1 : B Avalanche duration, T A v a l a n c h e du r a t i o nd i s t ., P ( T ) -4 -3 -2 -1 E = 1 : E = 1 :
185 mV E = 13 mV T P ( T ) -4 -2 L = 20 L = 40 L = 80 L = 99 E = 1 :
185 mV A L = 99 FIG. 6: Silent periods avalanche duration distributions (approach II). A: P ( T ) for E in threeregimes: ( (cid:72) ) subcritical, ( (cid:108) ) critical, and ( (cid:78) ) supercritical. Notice that the three regimes havePL shape inside the shaded area. A inset: critical regime P ( T ) for critical E = 1 .
185 mV anddifferent L . B: F ( T ) for critical E = 1 .
185 mV and different L used to fit Eq. (11) yielding T c ( L )and τ = 1 .
57. Solid line: example fit of Eq. (11) for L = 20 (R-Squared: 99 . T c ∼ exp( cL D ) yielding D = 1 . D of the avalanche size cutoff, s c (Fig. 5B inset), reveals the dimen-sionality of the path that propagates the avalanches [45]: D = 0 in the subcritical regimemeans that only small avalanches that remain confined in the bulk of the network (inde-pendently of L ) exist; D ≈ D ≈ D = 2 in the supercriticalregime means that avalanches are either spreading both radially and inside the columns oronly radially as a surface wave [45, 49].Figs. 5A and 6A show that there is indeed a PL shape for small avalanche sizes andduration for every considered E (either subcritical, critical or supercritical). The insetof these figures show that the cutoff of these distributions scale with system size L forcritical E = 1 .
185 mV. Figs. 5B and 6B show the cumulative distributions for the critical E = 1 .
185 mV and many L (the same cases as the insets of panels A in both figures), inwhich it is easier to notice the cutoff diverging as the system is increased. We fit Eq. (11) toevery cumulative distribution in Figs. 5B and 6B to estimate the cutoffs s c and T c . Thesecutoffs are plotted against L in the inset of both panels B of the respective figures. s c presentsthe usual behavior described by Eq. (7) with D = 0 . T c presents a stretchedexponential scaling with L , T c ∼ exp( kL D ) ( k is a constant) with D = 1 . . ≤ E ≤ . C. Power spectra
The power spectrum S ( f ) is related to temporal correlations via a Fourier transform.Thus, having a S ∼ /f means that there is long-range temporal correlations in the originalsignal. Such spectral density is found for the healthy brain state using Magnetoencephalog-raphy and Electroencephalography techniques [31, 33–36].The time series of avalanche sizes is constructed as follows: for approach I, the timeseries is generated by a sequence of stimuli, such that each stimulus generates a singleavalanche that spreads through a different realization of network disorder; for approach II,the sequence of avalanche sizes is given by the sequence of the sums of the activity between14wo subsequent silent intervals in Fig. 2. We calculated S ( f ) for the time series of avalanchesizes on the subcritical, critical and supercritical state. The results are shown in Fig. 7 forboth approaches of measuring avalanches.Absorbing state avalanches are generated by different stimuli on independent networkdisorder trials. This approach is then expected to lack the 1 /f profile because subsequentavalanches are statistically independent. Indeed, Fig. 7A shows that all the considered E have the exact same constant power spectrum, corresponding to white noise.On the other hand, silent periods avalanches are temporally dependent of one another,because they follow from the dynamics of the system. Fig. 7B shows the power spectrum forthis second approach averaged over many network disorder trials. It is clear that the powerspectrum follows S ∼ /f b with b ≈ . E = 1 .
185 mV. The subcriticalsystem also displays
S ∼ /f b , but has b (cid:38) . Frequency, f (Hz) P o w e r s p ec t r u m , S ( f ) E = 1 : E = 1 : E = 13mV1 =f b ; b = 1 : B Frequency, f -2 -1 P o w e r s p ec t r u m , S ( f ) E = 1 : E = 1 : E = 13mV A (Hz) FIG. 7: Power spectra of avalanche size time series for both approaches and three regimes of E : ( (cid:72) )subcritical, ( (cid:108) ) critical, and ( (cid:78) ) supercritical. A: S ( f ) for absorbing state avalanches (approachI). The three regimes have white noise power spectrum. B: S ( f ) for silent periods avalanches(approach II). The critical regime has S ∼ /f b with b ∼ . /f b with b = 0). In fact, within this approach, S ( f ) has been shown to display1 /f b with b varying from b (cid:38) . b = 0 while E varies from E = 1 . E = 13 mV [45]. IV. CONCLUDING REMARKS
We studied a model for the visual cortex of mammals consisting of a layered-columnarnetwork of compartmental neurons. The connections were randomly distributed accordingto structural parameters obtained from the literature [45]. This system presents an ab-sorbing phase transition from an inactive phase to a percolating phase with the adjustedparameters [45]. The spreading of the signal is deterministic and advances one neuronalcompartment per time step. Thus, apparent inactive intervals between peaks of activityemerge because we measure only soma spikes. This dynamics makes the avalanche defini-tion ambiguous: one may either define an avalanche as all the activity sparked by a singlestimulus between consecutive absorbing states or as all the activity between two silent pe-riods in the time series data. We acknowledge that the second approach is more similarto the experimental approach for measuring local field potentials avalanches. We use theEPSP parameter, E , in order to study the avalanche distributions in three different regimes:subcritical, critical and supercritical.Both approaches yield similar results in the critical region of the model: the avalanchesare PL distributed both in size and duration with similar exponents, contrary to whathas been found for another model [44]. Even though the critical behavior is similar, theavalanche distribution of both cases has different features with respect to the cutoff: whileapproach I yields a wider PL range for avalanche sizes and duration, the PL from approachII presents accumulation of large events on the cutoff bump [3]. This highlights the strongdependence of the results on the model details. The cutoff of these distributions definethe dimensionality of the avalanches as expected. However, the cutoffs of the distributionsof silent period avalanches scale with the system size similarly to the expected for criticalsystems even on the supercritical regime. Table II shows the values of all the exponents ofthe model for both approaches. The reported exponents have been verified via statisticaltests yielding no more than 7% error in their values. Details of the tests may be foundelsewhere [45]. 16 ABLE II: Summary of the model fitted exponents for the absorbing state avalanches (approachI) and for the silent period avalanches (approach II). Notice that the exponents defined here are notmutually exclusive because quantities have different interpretation in each approach. For instance,the duration cutoff in approach I is simply related to the dimensionality of avalanches (becausethe total spreading of the signal is constrained by the spatial extent of the network), whereasthe duration cutoff in approach II is the autocorrelation characteristic time [45]. Moreover, ¯ τ is related to the dynamical exponent inside the Griffiths phase and is non-universal [51]. Theexponents below have been verified through statistical tests and yielded no more than 7% error inthe reported values. Details of the test may be found elsewhere [45]. Quantity Approach I Approach II P ( s ) ∼ s − α Eq. (5) 1.4(1) 1.4(1) P ( T ) ∼ T − τ Eq. (6) 1.4(1) 1.5(1) s c ∼ L D Eq. (7) 1.7(1) 0.9(1)¯ T c ∼ L ¯ µ Eq. (8) 0.9(1) - T c ∼ exp( kL D ) - - 1.0(1) a = ( τ − / ( α −
1) - 1.1(1) 1.5(1)
S ∼ f − b - 0 1.2(1) Nevertheless, the results are significantly different in the subcritical and supercriticalregimes: the absorbing state avalanches do not present any PL distribution whereas the silentperiod avalanches presents a restricted PL shape for the distribution of small avalanches bothin size and duration (the power law does not scale with system size). This result is alsodependent on model details, given that other authors found no restricted power-laws outsideof the critical range using another kind of model [42, 43].The power spectra are also very different in the critical region of the model. This isexpected due to the strong dependence of the power spectrum on the avalanche definition.The independent avalanches generated by approach I results in a typical white noise powerspectrum, whereas approach II yields 1 /f b with b ≈ .
2. These results are in agreement withthe power spectra of other models that rely on the inhibitory-excitatory synaptic balancemechanism [24, 37]. 17bsorbing state avalanches are more precise for finding the critical region of the model.However, there are no guarantees that an experimental system will reach an absorbing statein between two measured avalanches. It is thus more likely that experimental PL avalanchedistributions that scale with system size alone will yield ambiguous results for defining theunderlying system critical behavior. Therefore, we emphasize that other quantities must beemployed in order to rigorously define the critical region of the system. For instance, wemay define an order parameter-susceptibility pair for this model and show that the orderparameter varies continuously to zero whilst the susceptibility diverges as the critical pointis approached [45]. Such procedure for finding the critical point is successfully applied in(non-)equilibrium statistical physics problems and is ultimately unambiguous.
Acknowledgments
We thank discussions with M. Copelli, R. Dickman and G. S. Bortolotto. MGS thankspartial support from CNPq during the development of this work, and funding from CIHRand the Mark Rayport and Shirley Ferguson Rayport Fellowship for Epilepsy Surgery of theMontreal Neurological Institute and Hospital. [1] P. Bak,
How Nature Works: The Science of Self-Organized Criticality (Copernicus Springer-Verlag New York, New York, NY, USA, 1996).[2] H. J. Jensen,
Self-Organized Criticality (Cambridge University Press, Cambridge, UK, 1998).[3] G. Pruessner,
Self-Organised Criticality (Cambridge University Press, Cambridge, UK, 2012).[4] V. Frette, K. Christensen, A. Malthe-Sørenssen, J. Feder, T. Jøssang, and P. Meakin, Nature , 49 (1996).[5] K. Christensen, H. Flyvbjerg, and Z. Olami, Phys. Rev. Lett. , 2737 (1993).[6] B. D. Malamud, G. Morein, and D. L. Turcotte, Science , 1840 (1998).[7] T. O. Richardson, E. J. H. Robinson, K. Christensen, H. J. Jensen, N. R. Franks, and A. B.Sendova-Franks, PLoS ONE , e9621 (2010), URL .[8] B. Gutenberg and C. F. Richter,
Seismicity Of The Earth And Associated Phenomena (Prince-ton University Press, Princeton, NJ, USA, 1949).
9] F. Y. Wang and Z. G. Dai, Nat. Phys. , 465 (2013), URL .[10] X. Gabaix, P. Gopikrishnan, V. Plerou, and H. E. Stanley, Nature , 267 (2003), URL .[11] R. Albert and A.-L. Barab´asi, Rev. Mod. Phys. , 47 (2002).[12] T. Gisiger, Biol. Rev. , 161 (2001), URL http://dx.doi.org/10.1017/S1464793101005607 .[13] D. R. Chialvo, Nat. Phys. , 744 (2010).[14] W. L. Shew and D. Plenz, Neuroscientist , 88 (2013), URL http://dx.doi.org/10.1177/1073858412445487 .[15] J. M. Beggs and N. Timme, Front. Physiol. , 163 (2012).[16] H. E. Stanley, Introduction to Phase Transitions and Critical Phenomena (Oxford UniversityPress, New York, UK, 1971).[17] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. , 381 (1987).[18] J. P. Sethna, K. A. Dahmen, and C. R. Myers, Nature , 242 (2001), URL .[19] H. J. Jensen, K. Christensen, and H. C. Fogedby, Phys. Rev. B , 7425 (1989).[20] D. Stassinopoulos and P. Bak, Phys. Rev. E , 5033 (1995).[21] M. Usher, M. Stemmler, and Z. Olami, Phys. Rev. Lett. , 326 (1995).[22] J. M. Beggs and D. Plenz, J. Neurosci. , 11167 (2003).[23] O. Kinouchi and M. Copelli, Nat. Phys. , 348 (2006).[24] L. de Arcangelis, C. Perrone-Capano, and H. J. Herrmann, Phys. Rev. Lett. , 028107 (2006).[25] A. Levina, J. M. Herrmann, and T. Geisel, Nat. Phys. , 857 (2007).[26] A. de Andrade Costa, M. Copelli, and O. Kinouchi, J. Stat. Mech. p. P06004 (2015), URL http://dx.doi.org/10.1088/1742-5468/2015/06/P06004 .[27] T. L. Ribeiro, M. Copelli, F. Caixeta, H. Belchior, D. R. Chialvo, M. A. L. Nicolelis, andS. Ribeiro, PLoS ONE , e14129 (2010).[28] M. Girardi-Schappo, O. Kinouchi, and M. H. R. Tragtenberg, Phys. Rev. E , 024701 (2013),URL http://dx.doi.org/10.1103/PhysRevE.88.024701 .[29] R. V. Williams-Garc´ıa, M. Moore, J. M. Beggs, and G. Ortiz, Phys. Rev. E , 062714(2014).[30] W. L. Shew, W. P. Clawson, J. Pobst, Y. Karimipanah, N. C. Wright, and R. Wessel, Nat. hys. , 659 (2015).[31] K. Linkenkaer-Hansen, V. V. Nikouline, J. M. Palva, and R. J. Ilmoniemi, J. Neurosci. ,1370 (2001).[32] A. Haimovici, E. Tagliazucchi, P. Balenzuela, and D. R. Chialvo, Phys. Rev. Lett. , 178101(2013).[33] E. Novikov, A. Novikov, D. Shannahoff-Khalsa, B. Schwartz, and J. Wright, Phys. Rev. E , R2387 (1997).[34] M. C. Teich, C. Heneghan, S. B. Lowen, T. Ozaki, and E. Kaplan, J. Opt. Soc. Am. A ,529 (1997).[35] J. Andrew Henrie and R. Shapley, J. Neurophysiol. , 479 (2005), URL http://dx.doi.org/10.1152/jn.00919.2004 .[36] D. Hermes, K. J. Miller, B. A. Wandell, and J. Winawer, Cereb. Cortex (2014).[37] F. Lombardi, H. J. Herrmann, and L. de Arcangelis, Chaos , 047402 (2017).[38] J. M. Beggs and D. Plenz, J. Neurosci. , 5216 (2004).[39] V. Priesemann, M. H. J. Munk, and M. Wibral, BMC Neurosci. , 40 (2009).[40] G. Hahn, T. Petermann, M. N. Havenith, S. Yu, W. Singer, D. Plenz, and D. Nikoli´c, Journalof Neurophysiology , 3312 (2010).[41] V. Priesemann, M. Wibral, M. Valderrama, R. Pr¨opper, M. Le Van Quyen, T. Geisel,J. Triesch, D. Nikoli´c, and M. H. J. Munk, Front Syst Neurosci. , 108 (2014), URL http://dx.doi.org/10.3389/fnsys.2014.00108 .[42] S. Scarpetta and A. de Candia, PLoS One , e64162 (2013), URL http://dx.doi.org/10.1371/journal.pone.0064162 .[43] S. Scarpetta and A. de Candia, Front. Syst. Neurosci. , 88 (2014), URL http://dx.doi.org/10.3389/fnsys.2014.00088 .[44] M. Martinello, J. Hidalgo, S. di Santo, A. Maritan, D. Plenz, and M. A. Mu˜noz,arXiv:1703.05079v1 [q-bio.NC] (2017).[45] M. Girardi-Schappo, G. S. Bortolotto, J. J. Gonsalves, L. T. Pinto, and M. H. R. Tragtenberg,Sci. Rep. , 29561 (2016), URL .[46] T. D. Albright, T. M. Jessell, E. R. Kandel, and M. I. Posner, Neuron , S1 (2000).[47] S. R. Williams and G. J. Stuart, Science , 1907 (2002).[48] M. E. J. Newman, Contemp. Phys. , 323 (2005).
49] G. S. Bortolotto, M. Girardi-Schappo, J. J. Gonsalves, L. T. Pinto, and M. H. R. Tragtenberg,J. Phys. Conf. Ser. , 012008 (2016), URL http://stacks.iop.org/1742-6596/686/i=1/a=012008 .[50] A. J. Bray and G. J. Rodgers, Phys. Rev. B , 11461 (1988).[51] R. Juh´asz, I. A. Kov´acs, and F. Igl´oi, Phys. Rev. E , 032815 (2015)., 032815 (2015).