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–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 , 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 , 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  or mechanisms (especially in living systems) other than criticality that can lead to power-law distributions [23–29]. In particular, Friedman and Landsberg  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  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)  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 . 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 .
Our analytical calculations are built on techniques originated from studies that shed light on structural properties of networks , spread of epidemics , properties of site percolation on undirected  and directed  networks, branching processes [32,53], spread of online information on Twitter , 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].
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
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].
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 
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
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 .
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 . It is worth mentioning that, for analyzing supercritical avalanches, methods based on Ref.  are also possible; nonetheless, such methods produce rather poor results .
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).
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 . 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 . 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 .
The authors thank J. M. Beggs for helpful comments on the manuscript. A.F. and J.P.G. acknowledge support from the Science Foundation Ireland (Grants No. 16/IA/4470 and No. 16/RC/3918). F.R. acknowledges support from the National Science Foundation (CMMI-1552487) and the U.S. Army Research Office (W911NF-16-1-0104).