key: cord-0010910-6w1ybmzh
authors: Faqeeh, Ali; Osat, Saeed; Radicchi, Filippo; Gleeson, James P.
title: Emergence of power laws in noncritical neuronal systems
date: 2019-07-02
journal: nan
DOI: 10.1103/physreve.100.010401
sha: 82781f926342d141427ab920f79d0f1bb374fff4
doc_id: 10910
cord_uid: 6w1ybmzh

Experimental and computational studies provide compelling evidence that neuronal systems are characterized by power-law distributions of neuronal avalanche sizes. This fact is interpreted as an indication that these systems are operating near criticality, and, in turn, typical properties of critical dynamical processes, such as optimal information transmission and stability, are attributed to neuronal systems. The purpose of this Rapid Communication is to show that the presence of power-law distributions for the size of neuronal avalanches is not a sufficient condition for the system to operate near criticality. Specifically, we consider a simplistic model of neuronal dynamics on networks and show that the degree distribution of the underlying neuronal network may trigger power-law distributions for neuronal avalanches even when the system is not in its critical regime. To certify and explain our findings we develop an analytical approach based on percolation theory and branching processes techniques.

In neuronal systems, but also in many other apparatuses crucial to living organisms, the emergence of power-law distributions [1] [2] [3] [4] [5] [6] has a remarkable importance. The unique form of distribution may indicate that the system is operating near a critical point [7] [8] [9] and is therefore benefitting from a series of potential advantages of critical systems [10] [11] [12] , such as optimum information transmission [1, 13] , dynamical range and sensitivity to sensory stimuli [14] , information capacity [15, 16] , and stability [13, 17] .

In neural systems, power-law distribution of avalanche sizes and durations (lifetimes) have been observed both in experiments [1, 4, 8] and computational models [17] [18] [19] . To validate that power-law distributions are indeed due to criticality, one needs to perform other tests [4] , including testing finite-size scaling relations [4, 20] and performing collapse of temporal profiles [4, 20, 21] . However, these techniques cannot determine what mechanisms or conditions keep or pose the system in the critical regime or if/how the system may lose its criticality. Nevertheless, other approaches can be used to demonstrate mechanisms or conditions that can lead to criticality of biological systems [22] or mechanisms (especially in living systems) other than criticality that can lead to power-law distributions [23] [24] [25] [26] [27] [28] [29] . In particular, Friedman and Landsberg [22] considered a simplistic model for neuronal dynamics and introduced a mechanism through which the hierarchical structure of neuronal networks can generate powerlaw distributions even far from criticality. The importance of network structure underlying neuronal dynamics for the * afaqeeh@iu.edu generation of power laws has been also reported in other studies [18, 19] .

In this Rapid Communication, we demonstrate that the degree distribution of the network underlying neural dynamics plays a fundamental role in the emergence of power-law distributions of avalanche sizes. To do so, we consider a simplified model of neural dynamics on networks, and show that, for some scale-free networks, avalanche sizes obey power-law distributions even in subcritical dynamical regime. Moreover, in other cases in which the avalanche size distribution is a power law with exponential cut-off, we disclose what structural parameters determine the cut-off size and show that even in such cases it is possible to observe distributions that are approximately power law over several orders of magnitude. In addition to numerical evidence, we provide an analytical description of the phenomenon relying on techniques borrowed from the theory of percolation [30] and branching processes [31] [32] [33] . We believe that our findings may have important implications in understanding properties of dynamics on realworld networks that have heavy-tailed degree distributions [2, 5, [34] [35] [36] [37] .

As mentioned above, we consider a simplistic model of neural avalanches for which we can show lucidly the impact of the network structure. In our model, an avalanche starts with a single activated neuron and, at each time step, every one of the active neurons fires a signal that stimulates all of their neighbors. This stimulus activates with a probability p each neighbor that has not been already activated. The avalanche of activities continues until no new neuron can be activated. This model is identical to the so-called independent cascade model, often considered in the context of opinion spreading in social networks [38] [39] [40] . For neural dynamics, it is a more realistic version of the Friedman-Landsberg model (FLM) [22] as in our model each time a neuron receives a stimulus it has the chance to become activated while in the FLM the activation does not depend on the number of stimulations. In spite of its simplicity, our model captures the fast timescale behavior of integrate-and-fire models [41, 42] . This fact follows from the simplifying assumptions that repetitive activation is neglected and the stimulations that activate a neuron (by increasing its potential to above its firing threshold) are set at random [19, 22, 43, 44] .

An advantage of this simplification is that our model is equivalent to a bond percolation model, thus, the avalanche size distribution is identical to the probability distribution π s that a randomly chosen node belongs to a percolation cluster of size s. This analogy enables us to consider a set of wellestablished techniques developed for percolation models and branching processes. In the following, we first describe our analytical calculations. Then, we show that our theoretical predictions are in very good agreement with the results of numerical simulations.

We provide a unifying framework that can describe the avalanche properties on both undirected and directed networks. We consider networks with negligible source-target correlation, i.e., the correlation between the degree values at the ending points of an edge. Nevertheless, for directed networks (DNs), we include analysis for networks with and without input-output correlation, i.e., the correlation between the values of indegree j and outdegree k of a node. Thus we consider three network types: undirected networks (UNs), uncorrelated directed network (UDNs), and input-output correlated directed networks (CDNs).

To generate UNs with specific degree distribution we use the configuration model [45] [46] [47] and for DNs we use an extended version of this model [48] . In particular, if the number of stubs of indegree and outdegree distributions are unbalanced, we remove, from a fraction u of nodes, some stubs of the distribution with more stubs such that its tail conserves its form [48] .

Our analytical calculations are built on techniques originated from studies that shed light on structural properties of networks [49] , spread of epidemics [50] , properties of site percolation on undirected [51] and directed [52] networks, branching processes [32, 53] , spread of online information on Twitter [33] , and relevant methods for obtaining the properties of generating functions [54, 55] . The findings most related to our calculations correspond to those of Refs. [51, 52] in which analytical results for the functional form of the distribution of cluster sizes in a site percolation process were reported. To improve upon the findings of these references, we substitute parts of the approaches they employed with our own techniques developed on the basis of Refs. [33, 49, 50, [53] [54] [55] .

For UNs, we consider the degree distribution p k of the network and for DNs we consider the indegree distribution P j , the outdegree distribution P k (note that we use k for the degree in UNs as well as for the outdegree in DNs as, we will later show that they play the same role in describing the avalanche sizes), and the joint degree distribution P jk which equals the fraction of nodes with indegree j and outdegree k. From the degree distributions we can obtain the excess degree distribution functions q k = k p k / k for UNs or q jk = jP jk / j for DNs which describe the probability that following a random edge we find a node with, respectively, degree k + 1 (UNs) or indegree j and outdegree k (DNs).

We are going to calculate the distribution of avalanche sizes π s . This quantity depends on ρ s , the probability that following an edge of the network we reach an avalanche (cluster) with size s [30, 49, 50, 56] . To calculate these quantities we will need to work with their generating functions defined as, respectively, H 0 (z) = ∞ s=1 π s z s and H 1 (z) = ∞ s=0 ρ s z s . We will also need the generating functions for degree k and the excess degree distributions, defined as

A UN on which a bond percolation process with occupation probability p is applied should be described by G n (z) = G n (1 − p + p z) instead [50, 57] , where n = 0, 1; it is straightforward to show that this property holds also for DNs. We use this fact to extend the governing equations that Newman derived for H 0 and H 1 in the absence of percolation [50] to our case; thus we get

The first difference between our calculations and the method of Refs. [51, 52] for calculation of cluster sizes is that we use the accurately derived Eqs. (3) and (4) instead of equations derived from heuristics [58] . The next steps of our approach include (i) calculation of the leading order nonanalytic behavior of H 0 (z) by finding the behavior of G 1 and G 0 around η = 1 − p + pH 1 (1), and (ii) using the asymptotic properties of generating functions [48, 54, 55 ] to obtain π s for large avalanche sizes (s 1) using the results of (i). To do so, we integrate the above equations with the methods described in [51, 52] and improve upon parts of these methods by combining them with techniques and ideas, including branching processes methods [33, 53] .

In these regimes, H 0 (1), which equals the probability that a randomly chosen node is in a finite cluster, can be set to 1 [for the supercritical regime, we can instead assume that H 0 (1) ≈ 1, if p is not much larger than the critical occupation probability p c ]; thus, according to Eq. (4) H 1 (1) = 1 too. Accordingly, to obtain the leading order behavior of H 1 (z) [from Eq. (3)] and H 0 (z) [from Eq. (4)], we can assume

1 and z . = 1 − w for w 1, and then expand the degree-dependent generating functions around z = 1 to get [48] 

as φ → 0, where p c is obtained analytically according to the results of Refs. [47, 51, 52] , and the terms with λ or λ are present only for scale-free networks; these effective exponents are

where λ, λ o , and λ i are the exponents for the tail of the distribution of, respectively, the degrees k in a UN, the outdegrees k in a DN, and the indegrees j in that DN. The other coefficients in Eqs. (5) and (6) depend on the network degree distribution [48] . Note that p c > 0 for the range of λ values we considered. We keep up to the third (second) leading order term of G 1 (G 0 ) and substitute the result in Eq. (3) [Eq. (4)] to get

and

where, in Eq. (7), δ = p − p c . Using Eqs. (7) and (8) we can show that the leading order nonanalytic term of H 0 (z), depending on the dynamical regime, has either the form Rw β or R √ 1 + s * w [48] , where β is a noninteger number and s * and R are constant. According to the asymptotic properties of generating functions [48, 54, 55] , the first form gives a π s with a power-law tail and the second form results in a power law with exponential decay. In particular, for λ > 3 and nonscale-free networks, π s = R 1 s −3/2 at the critical point (p c ) and π s = R 1 s −3/2 e −s/s * in the noncritical phases; however, for 2 < λ < 3, π s = R 2 s −[1+1/( λ−1)] at p c (see Supplemental Material [48] for the definitions obtained for R 1 and R 2 [59] ). For these cases, Refs. [51, 52] reported the same results for the functional form of π s ; however, the equations they used [instead of Eqs. (3) and (4)] underestimate the prefactors [48] . We also retrieve the result s * ∝ δ −2 calculated previously for UNs (using another method [51] ) and we discover that a similar relation also holds for DNs; furthermore, we find that the exponential decay factor s * is also controlled by the skewness of the degree distribution according to s * ≈ 2p 2 k k 3 δ 2 k 2 2 for UNs and s * ≈ 2p 2 k jk 2 δ 2 jk 2 for DNs. This indicates that even for skewed non-scale-free networks it is possible that a power-law distribution of avalanche sizes, expanded for several orders of magnitudes (i.e., as long as s s * ), emerges. On the other hand, despite the expectations of Refs. [51, 52] , at the subcritical regime of 2 < λ < 3, we get pure power-law distribution in the form

where R = a * (1 − u) , and a * . = 1/ k k −λ for k in the tail of the corresponding k distribution (degree or outdegree distribution) and u < 1 depends on the mass of that tail. To obtain this result we assumed that the solution of Eq. (7) for 2 < λ < 3 has the form φ ∼ a 1 w α 1 + a 2 w α 2 + · · · , where α 1 < α 2 < · · · , and used the dominant balance method [33, 44, 60] to obtain the correct form for the leading order terms [48] . As demonstrated by Eq. (9), the exponent of this distribution in UDNs is λ = λ = λ o , and in CDNs and UNs is λ; in CDNs the leading order term can be corrected with a λ order term if the difference between the exponents λ and λ is considerably small. Figure 1 shows that our predictions for the subcritical and critical regimes of 2 < λ < 3 capture very well the powerlaw behaviors of the tail of π s . In the procedure for deriving Eq. (9) we made no assumption about the sign of δ; however, we immediately notice that Eq. (9) is only valid for the subcritical regime since in the supercritical regime (that δ > 0) the prefactor and hence the probability π s will not be a real non-negative value. In with the same maximum s as that of the numerical data. The nearly vertical cut-offs in s are caused by the cumulative sum on data with finite maximum s and not by the shape of π s . In each panel, the numerical results for lower p values are those located at lower positions. The networks are constructed using our extension of the configuration model [48] . the next section, we show that a set of equations other than Eqs. (5) and (6) and a modified method should be used to obtain valid results for the supercritical regime of 2 < λ < 3.

We first consider the governing Eqs. (3) and (4) . Then as we know that in the supercritical phase H 1 (1) < 1, for z values close to 1 (or z . = 1 − w with w 1) we can write

where h, η < 1 and 1. Thus, in the right-hand side of Eqs. (3) and (4),

= z G 1 (η) + G 1 (η)p + 1 2 G 1 (η)p 2 2 + · · · . (14) Therefore, according to Eqs. (10) and (14),

where a 0 = G 1 (η) = h [see Eqs. (3) and (10)], a 1 = G 1 (η)p, and a 2 = 1 2 G 1 (η)p 2 . Equation (15) gives

Now we consider that the supercritical properties of H 1 (z) can be well approximated using the behavior of near the branch point z = 1; around this point we have

where s * = 4a 0 a 2 (1−a 1 ) 2 . Now,

where b = pG 0 (η) (1−a 1 ) 2a 2

. Then, according to the asymptotic properties of generating functions [48, 54, 55] ,

where s * . =

(1−pG 1 (η)) 2 and η is calculated using Eqs. (3) and (12) according to the prescription described in Sec. S3.2.d of [48] . Figures 2(e) and 2(f) show that Eq. (22) performs well in describing the distribution of avalanches with finite (nonextensive) sizes in the supercritical regime of 2 < λ < 3. It is worth noting that, in the supercritical regime, extremely large avalanches do also exist; such avalanches have a size that scales linearly with the network size. Hence, in the thermodynamic limit where the network size N → ∞, their Fig. 1(a) . At the supercritical phase, a bump appears at the tail of the distribution which corresponds to the percolating clusters whose sizes diverge at the thermodynamic limit (i.e., as N → ∞). (e), (f) The distribution s (the cumulative distribution of π s ) of the finite clusters (i.e., excluding the bump at the tail) for the supercritical regime of (e) a UN with λ = 3.3, p c = 0.23, and 5 × 10 7 nodes and (f) the UDN of Fig. 1(a) . The dots represent numerical simulations and the dashed lines are the theoretical results. Closer to the critical point a better agreement between theory and numerics is observed. In panels (e) and (f), the results for the lower p value are those located at a lower position. size also diverges. In finite networks, the effect of such avalanches on the distribution of avalanche sizes can be observed as a bump (in DNs) or a single point (in UNs) in the tail of the distribution. As the system moves further from the critical point this bump (or point in UNs) separates and moves away from the rest of the distribution [Figs. 2(a)-2(d)].

An interesting outcome of Eq. (22) is that the dependence of s * (the exponential decay parameter) on the inverse of δ is no longer purely quadratic; nonetheless, s * is still determined by p and skewness of degree distributions through a function that depends on p and the properties of G 1 at η [61] . It is worth mentioning that, for analyzing supercritical avalanches, methods based on Ref. [53] are also possible; nonetheless, such methods produce rather poor results [48] .

As we mentioned earlier, a prominent implication of our results for the noncritical cases is that even non-scale-free networks can produce avalanches distributed according to a power law for several orders of magnitude. This has significant implications for the experiments of neuronal dynamics that are commonly performed on small size samples [62, 63] ; this is because in such cases pure power laws and the powerlaw part of noncritical systems may be indistinguishable (see Fig. 3 ).

FIG. 3. The distribution π s for the subcritical regime of (a) CDNs and (b) UNs with 5000 nodes for 10 5 avalanches in each network. The curves, from top to bottom, correspond to a CDN (UN) with, respectively, exponential, power-law, and log-normal outdegree (degree) distribution. For each of these networks an apparent powerlaw tail is observed for π s . (Note that power laws observed in experimental setups have similar ranges of s [1, 8] .) The simulations are performed far from the critical point at occupation probabilities in the range (0, p c /2). In non-scale-free CDNs, the correlations are implemented by setting j = k 0.7 .

In summary, we have provided analytical proofs, accompanied by numerical confirmations, that even in a simplified description of neuronal dynamics, because of structural heterogeneity, different types of critical and noncritical powerlaw avalanches with different exponents can be observed. This finding may help us to better explain the emergence of power-law neuronal avalanches with exponents different from 3/2 observed in experiments [64, 65] and in realistic computational simulations [19] . Moreover, critical systems are known to have crucial advantages such as optimum information transmission, capacity, and stability [1, 13, [15] [16] [17] because of their power-law avalanches. Thus the existence of power laws in other dynamical regimes implies that some noncritical systems may benefit from similar advantages due to divergence of the mean values and scale invariance of their power-law distributions [13] . Furthermore, the emergence of such noncritical power laws introduces new challenges for accurate detection of criticality in experimental setups due to the finite size of the commonly used samples. In addition to the importance of these findings in the context of neuronal systems, the insights on percolation properties of networks that this Rapid Communication provides may find applications in topics such as network robustness [56, 66] , epidemic spreading [50, 67, 68] , and stability of biological systems [69] .

Neuronal avalanches in neocortical circuits

Ising-like dynamics in large-scale functional brain networks

Self-organized criticality in living systems

Universal Critical Dynamics in High Resolution Neuronal Avalanche Data

Power-law distributions in empirical data

Colloquium: Criticality and dynamical scaling in living systems

Emergent complex neural dynamics

Being critical of criticality in the brain

Criticality in Neural Systems

Computation at the edge of chaos: Phase transitions and emergent computation

How Nature Works: The Science of Self-Organized Criticality

Critical brain networks

The criticality hypothesis: How local cortical networks might optimize information processing

Optimal dynamical range of excitable networks at criticality

Information capacity and transmission are maximized in balanced cortical networks with neuronal avalanches

Real-time computation at the edge of chaos in recurrent neural networks

Critical Branching Captures Activity in Living Neural Networks and Maximizes the Number of Metastable States

Recurrent infomax generates cell assemblies, neuronal avalanches, and simple cell-like selectivity

Neurobiologically realistic determinants of self-organized criticality in networks of spiking neurons

Crackling noise

Temporal profiles of avalanches on networks

Hierarchical networks, power laws, and neuronal avalanches

Power laws, Pareto distributions, and Zipf's law

Critical truths about power laws

Neutral Theory and Scale-Free Neural Dynamics

Avalanches in a stochastic model of spiking neurons

Power-law statistics and universal scaling in the absence of criticality

Can power-law scaling and neuronal avalanches arise from stochastic dynamics?

Can a time varying external drive give rise to apparent criticality in neural systems?

Networks: An Introduction

Branching Processes

Sandpile on Scale-Free Networks

Competition-Induced Criticality in a Model of Meme Popularity

Evolution of the social network of scientific collaborations

Functional clusters, hubs, and communities in the cortical microconnectome

Scale-free networks: A decade and beyond

Scale-free networks are rare

Talk of the network: A complex systems look at the underlying process of word-ofmouth

Using complex systems analysis to advance marketing theory development: Modeling heterogeneity effects on new product growth through stochastic cellular automata

Maximizing the spread of influence through a social network

A review of the integrate-and-fire neuron model: I. Homogeneous synaptic input

A review of the integrate-and-fire neuron model: II. Inhomogeneous synaptic input and network properties

Dynamical synapses causing self-organized criticality in neural networks

Percolation and its relations to other processes in networks

The asymptotic number of labeled graphs with given degree sequences

Complex networks: Structure and dynamics

A critical point for random graphs with a given degree sequence

for more details of calculations and/or further relevant results

Random graphs with arbitrary degree distributions and their applications

Spread of epidemic disease on networks

Percolation critical exponents in scale-free networks

Percolation in directed scale-free networks

Scale-free random branching trees in supercritical phase

Generating Functionology (A K Peters

Aspects and Applications of the Random Walk

Network Robustness and Fragility: Percolation on Random Graphs

Component sizes in networks with arbitrary degree distributions

G n (1 − p + p + p z) (that holds for site percolation too) is used in Refs

Note that obtaining the prefactors is important as the full functional form is needed to reliably compare the analytical predictions with the simulations results; moreover

Applied Asymptotic Analysis, Graduate studies in mathematics

Note that in other cases s * depended on the properties of G 1 at 1 which we could describe them in terms of averages of the different powers of k (and also j in DNs)

The use of neuronal networks on multielectrode arrays as biosensors

Highresolution CMOS MEA platform to study neurons at subcellular, cellular, and network levels

Spike avalanches exhibit universal dynamics across the sleep-wake cycle

Criticality in large-scale brain FMRI dynamics unveiled by a novel point process analysis

Emergence of coexisting percolating clusters in networks

Network cloning unfolds the effect of clustering on dynamical processes

Griffiths Phases on Complex Networks

Dynamical Instability in Boolean Networks as a Percolation Problem