Emergence of power laws in noncritical neuronal systems.

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.

INTRODUCTION

In neuronal systems, but also in many other apparatuses crucial to living organisms, the emergence of power-law distributions [1-6] has a remarkable importance. The unique form of distribution may indicate that the system is operating near a critical point [7-9] and is therefore benefitting from a series of potential advantages of critical systems [10-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-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-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 power-law distributions even far from criticality. The importance of network structure underlying neuronal dynamics for the 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-33]. We believe that our findings may have important implications in understanding properties of dynamics on real-world networks that have heavy-tailed degree distributions [2,5,34-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 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-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 that a randomly chosen node belongs to a percolation cluster of size . This analogy enables us to consider a set of well-established 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 and outdegree 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-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 of nodes, some stubs of the distribution with more stubs such that its tail conserves its form [48].

RESULTS

Relevant theoretical findings

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-55].

Constructing the governing equations

For UNs, we consider the degree distribution of the network and for DNs we consider the indegree distribution , the outdegree distribution (note that we use 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 which equals the fraction of nodes with indegree and outdegree . From the degree distributions we can obtain the excess degree distribution functions for UNs or for DNs which describe the probability that following a random edge we find a node with, respectively, degree (UNs) or indegree and outdegree (DNs).

We are going to calculate the distribution of avalanche sizes . This quantity depends on , the probability that following an edge of the network we reach an avalanche (cluster) with size [30,49,50,56]. To calculate these quantities we will need to work with their generating functions defined as, respectively, and . We will also need the generating functions for degree and the excess degree distributions, defined as A UN on which a bond percolation process with occupation probability is applied should be described by instead [50,57], where ; 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 and 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 by finding the behavior of and around , and (ii) using the asymptotic properties of generating functions [48,54,55] to obtain for large avalanche sizes () 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].

Solution methods for different regimes of dynamics

The critical and subcritical regimes

In these regimes, , 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 , if is not much larger than the critical occupation probability ]; thus, according to Eq. (4) too. Accordingly, to obtain the leading order behavior of [from Eq. (3)] and [from Eq. (4)], we can assume , where and for , and then expand the degree-dependent generating functions around to get [48] as , where 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 , , and are the exponents for the tail of the distribution of, respectively, the degrees in a UN, the outdegrees in a DN, and the indegrees in that DN. The other coefficients in Eqs. (5) and (6) depend on the network degree distribution [48]. Note that for the range of values we considered.

We keep up to the third (second) leading order term of () and substitute the result in Eq. (3) [Eq. (4)] to get and where, in Eq. (7), . Using Eqs. (7) and (8) we can show that the leading order nonanalytic term of , depending on the dynamical regime, has either the form or [48], where is a noninteger number and and are constant. According to the asymptotic properties of generating functions [48,54,55], the first form gives a with a power-law tail and the second form results in a power law with exponential decay. In particular, for and non-scale-free networks, at the critical point () and in the noncritical phases; however, for , at (see Supplemental Material [48] for the definitions obtained for and [59]). For these cases, Refs. [51,52] reported the same results for the functional form of ; however, the equations they used [instead of Eqs. (3) and (4)] underestimate the prefactors [48]. We also retrieve the result 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 is also controlled by the skewness of the degree distribution according to for UNs and 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 ), emerges.

On the other hand, despite the expectations of Refs. [51,52], at the subcritical regime of , we get pure power-law distribution in the form where , and for in the tail of the corresponding distribution (degree or outdegree distribution) and depends on the mass of that tail. To obtain this result we assumed that the solution of Eq. (7) for has the form , where , 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 , 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 capture very well the power-law behaviors of the tail of .

FIG. 1.

The probability and its cumulative distribution for a UDN with and (top panels) and a CDN with and (bottom panels) both with nodes. The green or blue dots are the numerical results and the dashed lines correspond to the theoretical power law [Eq. (9) for ] with the same maximum as that of the numerical data. The nearly vertical cut-offs in are caused by the cumulative sum on data with finite maximum and not by the shape of . In each panel, the numerical results for lower values are those located at lower positions. The networks are constructed using our extension of the configuration model [48].

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 ) the prefactor and hence the probability will not be a real non-negative value. In 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 .

The supercritical regime of

We first consider the governing Eqs. (3) and (4). Then as we know that in the supercritical phase , for values close to 1 (or with ) we can write where and . Thus, in the right-hand side of Eqs. (3) and (4), where . Now, we can write Therefore, according to Eqs. (10) and (14), where [see Eqs. (3) and (10)], , and . Equation (15) gives Now we consider that the supercritical properties of can be well approximated using the behavior of near the branch point ; around this point we have where . Now, where . Then, according to the asymptotic properties of generating functions [48,54,55], where 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 . 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 , their 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)].

FIG. 2.

(a)–(d) The full cluster size distribution at the critical and supercritical regimes of our neuronal dynamics model for a UDN with 5000 nodes and the other parameters identical to those of 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 ). (e), (f) The distribution (the cumulative distribution of ) of the finite clusters (i.e., excluding the bump at the tail) for the supercritical regime of (e) a UN with , , and 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 value are those located at a lower position.

An interesting outcome of Eq. (22) is that the dependence of (the exponential decay parameter) on the inverse of is no longer purely quadratic; nonetheless, is still determined by and skewness of degree distributions through a function that depends on and the properties of 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 power-law part of noncritical systems may be indistinguishable (see Fig. 3).

FIG. 3.

The distribution for the subcritical regime of (a) CDNs and (b) UNs with 5000 nodes for 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 power-law tail is observed for . (Note that power laws observed in experimental setups have similar ranges of [1,8].) The simulations are performed far from the critical point at occupation probabilities in the range . In non-scale-free CDNs, the correlations are implemented by setting .

CONCLUSIONS

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 power-law 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-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].