Physical Review. E
American Physical Society
Emergence of power laws in noncritical neuronal systems
Volume: 100, Issue: 1
DOI 10.1103/PhysRevE.100.010401
•
•
•
• Altmetric

### Notes

Abstract

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.

Faqeeh, Osat, Radicchi, and Gleeson: Emergence of power laws in noncritical neuronal systems

## I.. 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 $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–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 ${\pi }_{s}$ that a randomly chosen node belongs to a percolation cluster of size $s$. 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 $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–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].

## II.. RESULTS

### A.. 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].

### B.. Constructing the governing equations

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}=j{P}_{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 ${\pi }_{s}$. This quantity depends on ${\rho }_{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}\left(z\right)={\sum }_{s=1}^{\infty }{\pi }_{s}{z}^{s}$ and ${H}_{1}\left(z\right)={\sum }_{s=0}^{\infty }{\rho }_{s}{z}^{s}$. We will also need the generating functions for degree $k$ and the excess degree distributions, defined as

$\begin{array}{c}\hfill {G}_{0}\left(z\right)=\left\{\begin{array}{cc}\hfill \phantom{\rule{0.16em}{0ex}}\sum _{k=0}^{\infty }{p}_{k}{z}^{k},& \hfill \phantom{\rule{0.16em}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\text{UN}\\ \\ \hfill \sum _{j,k=0}^{\infty }{P}_{jk}{z}^{k},& \hfill \phantom{\rule{0.16em}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\text{DN},\end{array}\right\\end{array}$
$\begin{array}{c}\hfill {G}_{1}\left(z\right)=\left\{\begin{array}{cc}\hfill \phantom{\rule{0.16em}{0ex}}\sum _{k=0}^{\infty }{q}_{k}{z}^{k},& \hfill \phantom{\rule{0.16em}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\text{UN}\\ \\ \hfill \sum _{j,k=0}^{\infty }{q}_{jk}{z}^{k},& \hfill \phantom{\rule{0.16em}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\text{DN}.\end{array}\right\\end{array}$
A UN on which a bond percolation process with occupation probability $p$ is applied should be described by ${\stackrel{^}{G}}_{n}\left(z\right)={G}_{n}\left(1-p+p\phantom{\rule{0.16em}{0ex}}z\right)$ 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
$\begin{array}{c}\hfill {H}_{1}\left(z\right)=z{G}_{1}\left[1-p+p{H}_{1}\left(z\right)\right],\end{array}$
$\begin{array}{c}\hfill {H}_{0}\left(z\right)=z{G}_{0}\left[1-p+p{H}_{1}\left(z\right)\right].\end{array}$
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}\left(z\right)$ by finding the behavior of ${G}_{1}$ and ${G}_{0}$ around $\overline{\eta }=1-p+p{H}_{1}\left(1\right)$, and (ii) using the asymptotic properties of generating functions [48,54,55] to obtain ${\pi }_{s}$ for large avalanche sizes ($s\gg 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].

### C.. Solution methods for different regimes of dynamics

#### 1.. The critical and subcritical regimes

In these regimes, ${H}_{0}\left(1\right)$, 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}\left(1\right)\approx 1$, if $p$ is not much larger than the critical occupation probability ${p}_{\mathrm{c}}$]; thus, according to Eq. (4)${H}_{1}\left(1\right)=1$ too. Accordingly, to obtain the leading order behavior of ${H}_{1}\left(z\right)$ [from Eq. (3)] and ${H}_{0}\left(z\right)$ [from Eq. (4)], we can assume ${H}_{1}\left(1-w\right)\sim 1-\varphi$, where $\varphi \ll 1$ and $z\doteq 1-w$ for $w\ll 1$, and then expand the degree-dependent generating functions around $z=1$ to get [48]

$\begin{array}{cc}& {G}_{1}\left(1-p\varphi \right)\sim {G}_{1}\left(1\right)-{G}_{1}^{\prime }\left(1\right)p\varphi +{G}_{1}^{″}\left(1\right){\left(p\varphi \right)}^{2}+D{\varphi }^{\stackrel{~}{\lambda }-1}\hfill \\ & \phantom{\rule{1em}{0ex}}+o\left({\varphi }^{2},{\varphi }^{\stackrel{~}{\lambda }-1}\right)=1-\frac{p}{{p}_{\mathrm{c}}}\varphi +B{\varphi }^{2}+D{\varphi }^{\stackrel{~}{\lambda }-1}\hfill \\ & \phantom{\rule{83.5pt}{0ex}}+o\left({\varphi }^{2},{\varphi }^{\stackrel{~}{\lambda }-1}\right),\hfill \end{array}$
$\begin{array}{cc}& {G}_{0}\left(1-p\varphi \right)\sim {G}_{0}\left(1\right)-{G}_{0}^{\prime }\left(1\right)p\varphi +M{\varphi }^{\overline{\lambda }-1}+\mathcal{O}\left({\varphi }^{2}\right)\hfill \\ & \phantom{\rule{1em}{0ex}}=\phantom{\rule{0.16em}{0ex}}1-E\varphi +M{\varphi }^{\overline{\lambda }-1}+\mathcal{O}\left({\varphi }^{2}\right),\hfill \end{array}$
as $\varphi \to 0$, where ${p}_{\mathrm{c}}$ is obtained analytically according to the results of Refs. [47,51,52], and the terms with $\overline{\lambda }$ or $\stackrel{~}{\lambda }$ are present only for scale-free networks; these effective exponents are
$\begin{array}{c}\hfill \overline{\lambda }=\left\{\begin{array}{cc}\hfill \lambda & \hfill \text{(UN)}\\ \hfill {\lambda }_{\mathrm{o}}& \hfill \text{(DN)}\end{array}\right\\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}\stackrel{~}{\lambda }=\left\{\begin{array}{cc}\hfill \lambda -1& \hfill \text{(UN)}\\ \hfill {\lambda }_{\mathrm{o}}& \hfill \text{(UDN)}\\ \hfill {\lambda }_{\mathrm{o}}-\frac{{\lambda }_{\mathrm{o}}-1}{{\lambda }_{\mathrm{i}}-1}& \hfill \text{(CDN)},\end{array}\right\\end{array}$
where $\lambda$, ${\lambda }_{\mathrm{o}}$, and ${\lambda }_{\mathrm{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}_{\mathrm{c}}>0$ for the range of $\stackrel{~}{\lambda }$ 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

$\begin{array}{cc}\hfill {H}_{1}\left(1-w\right)& \doteq 1-\varphi \sim \left(1-w\right)\left(1-\frac{p}{{p}_{\mathrm{c}}}\varphi +B{\varphi }^{2}+D{\varphi }^{\stackrel{~}{\lambda }-1}\right)\hfill \\ & ⇒\phantom{\rule{0.16em}{0ex}}\phantom{\rule{4pt}{0ex}}w\sim -\frac{\delta }{{p}_{\mathrm{c}}}\varphi +B{\varphi }^{2}+D{\varphi }^{\stackrel{~}{\lambda }-1}\hfill \end{array}$
and
$\begin{array}{cc}\hfill {H}_{0}\left(1-w\right)\sim & \phantom{\rule{0.16em}{0ex}}1-E\varphi +M{\varphi }^{\overline{\lambda }-1},\hfill \end{array}$
where, in Eq. (7), $\delta =p-{p}_{\mathrm{c}}$. Using Eqs. (7) and (8) we can show that the leading order nonanalytic term of ${H}_{0}\left(z\right)$, depending on the dynamical regime, has either the form $R{w}^{\phantom{\rule{0.16em}{0ex}}\beta }$ or $R\sqrt{1+{s}^{*}w}$[48], where $\beta$ 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 ${\pi }_{s}$ with a power-law tail and the second form results in a power law with exponential decay. In particular, for $\stackrel{~}{\lambda }>3$ and non-scale-free networks, ${\pi }_{s}={R}_{1}{s}^{-3/2}$ at the critical point (${p}_{\mathrm{c}}$) and ${\pi }_{s}={R}_{1}{s}^{-3/2}{e}^{-s/{s}^{*}}$ in the noncritical phases; however, for $2<\stackrel{~}{\lambda }<3$, ${\pi }_{s}={R}_{2}{s}^{-\left[1+1/\left(\stackrel{~}{\lambda }-1\right)\right]}$ at ${p}_{\mathrm{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 ${\pi }_{s}$; however, the equations they used [instead of Eqs. (3) and (4)] underestimate the prefactors [48]. We also retrieve the result ${s}^{*}\propto {\delta }^{-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}^{*}\approx \frac{2{p}^{2}〈k〉〈{k}^{3}〉}{{\delta }^{2}{〈{k}^{2}〉}^{2}}$ for UNs and ${s}^{*}\approx \frac{2{p}^{2}〈k〉〈j{k}^{2}〉}{{\delta }^{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\ll {s}^{*}$), emerges.

On the other hand, despite the expectations of Refs. [51,52], at the subcritical regime of $2<\stackrel{~}{\lambda }<3$, we get pure power-law distribution in the form

$\begin{array}{c}\hfill {\pi }_{s}\sim \left\{\begin{array}{cc}\hfill \stackrel{^}{R}\left(1+\frac{p{p}_{\mathrm{c}}〈k〉}{-\delta }\right){\left(\frac{p{p}_{\mathrm{c}}}{-\delta }\right)}^{{\lambda }_{\mathrm{o}}-1}{s}^{-{\lambda }_{\mathrm{o}}}& \hfill \text{(UDN)}\\ \\ \hfill \stackrel{^}{R}{\left(\frac{p{p}_{\mathrm{c}}}{-\delta }\right)}^{\overline{\lambda }-1}\left[{\left(\frac{p{p}_{\mathrm{c}}}{-\delta }\right)}^{\stackrel{~}{\lambda }-\overline{\lambda }+1}{s}^{-\stackrel{~}{\lambda }}+{s}^{-\overline{\lambda }}\right]& \hfill \text{(CDN/UN)},\end{array}\right\\end{array}$
where $\stackrel{^}{R}={a}^{*}\left(1-\overline{u}\right)$, and ${a}^{*}\doteq 1/{\sum }_{\overline{k}}{\overline{k}}^{\phantom{\rule{0.16em}{0ex}}-\overline{\lambda }}$ for $\overline{k}$ in the tail of the corresponding $k$ distribution (degree or outdegree distribution) and $\overline{u}<1$ depends on the mass of that tail. To obtain this result we assumed that the solution of Eq. (7) for $2<\stackrel{~}{\lambda }<3$ has the form $\varphi \sim {a}_{1}{w}^{{\alpha }_{1}}+{a}_{2}{w}^{{\alpha }_{2}}+\cdots$, where ${\alpha }_{1}<{\alpha }_{2}<\cdots$, 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 $\stackrel{~}{\lambda }=\overline{\lambda }={\lambda }_{\mathrm{o}}$, and in CDNs and UNs is $\stackrel{~}{\lambda }$; in CDNs the leading order term can be corrected with a $\overline{\lambda }$ order term if the difference between the exponents $\stackrel{~}{\lambda }$ and $\overline{\lambda }$ is considerably small. Figure 1 shows that our predictions for the subcritical and critical regimes of $2<\stackrel{~}{\lambda }<3$ capture very well the power-law behaviors of the tail of ${\pi }_{s}$.
FIG. 1.
The probability ${\pi }_{s}$ and its cumulative distribution ${\mathrm{\Omega }}_{s}$ for a UDN with ${\lambda }_{\mathrm{o}}=2.7$ and ${\lambda }_{\mathrm{i}}=3.7$ (top panels) and a CDN with ${\lambda }_{\mathrm{o}}=3.1$ and ${\lambda }_{\mathrm{i}}=3.7$ (bottom panels) both with $2×{10}^{7}$ nodes. The green or blue dots are the numerical results and the dashed lines correspond to the theoretical power law ${\pi }_{s}$ [Eq. (9) for $p<{p}_{\mathrm{c}}$] with the same maximum $s$ as that of the numerical data. The nearly vertical cut-offs in ${\mathrm{\Omega }}_{s}$ are caused by the cumulative sum on data with finite maximum $s$ and not by the shape of ${\pi }_{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].

In the procedure for deriving Eq. (9) we made no assumption about the sign of $\delta$; however, we immediately notice that Eq. (9) is only valid for the subcritical regime since in the supercritical regime (that $\delta >0$) the prefactor and hence the probability ${\pi }_{s}$ 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 $2<\stackrel{~}{\lambda }<3$.

#### 2.. 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}\left(1\right)<1$, for $z$ values close to 1 (or $z\doteq 1-w$ with $w\ll 1$) we can write

${H}_{1}\left(z\right)=\stackrel{~}{h}+\epsilon =1-\eta +\epsilon ,$
where $\stackrel{~}{h},\eta <1$ and $\epsilon \ll 1$. Thus, in the right-hand side of Eqs. (3) and (4),
$1-p+p{H}_{1}\left(z\right)=1-p\eta +p\epsilon$
$\doteq \overline{\eta }+p\epsilon ,$
where $\overline{\eta }=1-p\eta =1-p\left(1-\stackrel{~}{h}\right)=1-p+p\stackrel{~}{h}$. Now, we can write
${H}_{1}\left(z\right)=z\phantom{\rule{0.16em}{0ex}}{G}_{1}\left(\overline{\eta }+p\epsilon \right)$
$=z\left[{G}_{1}\left(\overline{\eta }\right)\phantom{\rule{0.16em}{0ex}}+\phantom{\rule{0.16em}{0ex}}{G}_{1}^{\prime }\left(\overline{\eta }\right)p\epsilon \phantom{\rule{0.16em}{0ex}}+\phantom{\rule{0.16em}{0ex}}\frac{1}{2}{G}_{1}^{″}\left(\overline{\eta }\right){p}^{2}{\epsilon }^{2}\phantom{\rule{0.16em}{0ex}}+\cdots \right].\phantom{\rule{2.em}{0ex}}$
Therefore, according to Eqs. (10) and (14),
$\stackrel{~}{h}+\epsilon =z\phantom{\rule{0.16em}{0ex}}\left[{a}_{0}+{a}_{1}\epsilon +{a}_{2}{\epsilon }^{2}+\cdots \right],$
where ${a}_{0}={G}_{1}\left(\overline{\eta }\right)=\stackrel{~}{h}$ [see Eqs. (3) and (10)], ${a}_{1}={G}_{1}^{\prime }\left(\overline{\eta }\right)p$, and ${a}_{2}=\frac{1}{2}{G}_{1}^{″}\left(\overline{\eta }\right){p}^{2}$. Equation (15) gives
$\epsilon \sim \frac{1-{a}_{1}z}{2{a}_{2}z}±\frac{1}{2{a}_{2}z}\sqrt{{\left(1-{a}_{1}z\right)}^{2}+4{a}_{0}{a}_{2}z\left(1-z\right)}.$
Now we consider that the supercritical properties of ${H}_{1}\left(z\right)$ can be well approximated using the behavior of $\epsilon$ near the branch point $z=1$; around this point we have
$\epsilon \sim \frac{1-{a}_{1}}{2{a}_{2}}±\frac{1-{a}_{1}}{2{a}_{2}}\sqrt{1+{s}^{*}\left(1-z\right)},$
where ${s}^{*}=\frac{4{a}_{0}{a}_{2}}{{\left(1-{a}_{1}\right)}^{2}}$. Now,
${H}_{0}\left(z\right)=z\phantom{\rule{0.16em}{0ex}}{G}_{0}\left(\overline{\eta }+p\epsilon \right)$
$\sim \left(1-w\right)\left[{G}_{0}\left(\overline{\eta }\right)+{G}_{0}^{\prime }\left(\overline{\eta }\right)p\epsilon +o\left({\epsilon }^{2}\right)\right]$
$\sim \text{analytical}\phantom{\rule{4.pt}{0ex}}\text{terms}±b\sqrt{1+{s}^{*}\left(1-z\right)},$
where $b=\frac{p{G}_{0}^{\prime }\left(\overline{\eta }\right)\phantom{\rule{0.16em}{0ex}}\left(1-{a}_{1}\right)}{2{a}_{2}}$. Then, according to the asymptotic properties of generating functions [48,54,55],
${\pi }_{s}\sim ±\frac{b{{s}^{*}}^{1/2}}{2\sqrt{\pi }}\phantom{\rule{0.222222em}{0ex}}{s}^{-3/2}\phantom{\rule{0.222222em}{0ex}}{e}^{-s/{s}^{*}}\phantom{\rule{0.16em}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\text{as}\phantom{\rule{0.16em}{0ex}}s\to \infty$
$\sim ±\frac{\phantom{\rule{0.16em}{0ex}}{G}_{0}^{\prime }\left(\overline{\eta }\right)\sqrt{{G}_{1}\left(\overline{\eta }\right)}}{\sqrt{2\pi \phantom{\rule{0.16em}{0ex}}{G}_{1}^{″}\left(\overline{\eta }\right)}}\phantom{\rule{0.222222em}{0ex}}{s}^{-3/2}\phantom{\rule{0.222222em}{0ex}}{e}^{-s/{s}^{*}}\phantom{\rule{0.16em}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\text{as}\phantom{\rule{0.16em}{0ex}}s\to \infty ,$
where ${s}^{*}\doteq \frac{2{p}^{2}{G}_{1}\left(\overline{\eta }\right)\phantom{\rule{0.16em}{0ex}}{G}_{1}^{″}\left(\overline{\eta }\right)}{{\left(1-p{G}_{1}^{\prime }\left(\overline{\eta }\right)\right)}^{2}}$ and $\overline{\eta }$ 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<\stackrel{~}{\lambda }<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\to \infty$, 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 $N\to \infty$). (e), (f) The distribution ${\mathrm{\Omega }}_{s}$ (the cumulative distribution of ${\pi }_{s}$) of the finite clusters (i.e., excluding the bump at the tail) for the supercritical regime of (e) a UN with $\lambda =3.3$, ${p}_{\mathrm{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.

An interesting outcome of Eq. (22) is that the dependence of ${s}^{*}$ (the exponential decay parameter) on the inverse of $\delta$ 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 $\overline{\eta }$[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 ${\pi }_{s}$ for the subcritical regime of (a) CDNs and (b) UNs with 5000 nodes for ${10}^{\phantom{\rule{0.16em}{0ex}}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 power-law tail is observed for ${\pi }_{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 $\left(0,{p}_{\mathrm{c}}/2\right)$. In non-scale-free CDNs, the correlations are implemented by setting $j={k}^{\phantom{\rule{0.16em}{0ex}}0.7}$.

## III.. 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].

## ACKNOWLEDGMENTS

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).

## REFERENCES

[1]

J. M. Beggs and D. Plenz, . Neuronal avalanches in neocortical circuits, J. Neurosci.23, 11167 (2003).JNRSDS0270-6474

[2]

D. Fraiman, P. Balenzuela, J. Foss, and D. R. Chialvo, . Ising-like dynamics in large-scale functional brain networks, Phys. Rev. E79, 061922 (2009).PLEEE81539-3755

[3]

C. Adami, . Self-organized criticality in living systems, Phys. Lett. A203, 29 (1995)., doi: 10.1016/0375-9601(95)00372-A

[4]

N. Friedman, S. Ito, B. A. W. Brinkman, M. Shimono, R. E. Lee DeVille, K. A. Dahmen, J. M. Beggs, and T. C. Butler, . Universal Critical Dynamics in High Resolution Neuronal Avalanche Data, Phys. Rev. Lett.108, 208102 (2012).PRLTAO0031-9007

[5]

A. Clauset, C. Shalizi, and M. Newman, . Power-law distributions in empirical data, SIAM Rev.51, 661 (2009).SIREAD0036-1445

[6]

M. A. Muñoz, . Colloquium: Criticality and dynamical scaling in living systems, Rev. Mod. Phys.90, 031001 (2018).RMPHAT0034-6861

[7]

D. R. Chialvo, . Emergent complex neural dynamics, Nat. Phys.6, 744 (2010).1745-2473, doi: 10.1038/nphys1803

[8]

J. M. Beggs and N. Timme, . Being critical of criticality in the brain, Front. Physiol.3, 163 (2012).1664-042X, doi: 10.3389/fphys.2012.00163

[9]

D. Plenz and E. Niebur, Criticality in Neural Systems (Wiley-VCH, Weinheim, Germany, 2014).

[10]

C. G. Langton, . Computation at the edge of chaos: Phase transitions and emergent computation, Phys. D (Amsterdam, Neth.)42, 12 (1990).PDNPDT0167-2789

[11]

P. Bak, How Nature Works: The Science of Self-Organized Criticality (Springer & Business Media, New York, 1996).

[12]

D. R. Chialvo, . Critical brain networks, Phys. A (Amsterdam, Neth.)340, 756 (2004).PHYADX0378-4371

[13]

J. M. Beggs, . The criticality hypothesis: How local cortical networks might optimize information processing, Philos. Trans. R. Soc. A366, 329 (2008).PTRMAD1364-503X

[14]

O. Kinouchi and M. Copelli, . Optimal dynamical range of excitable networks at criticality, Nat. Phys.2, 348 (2006).1745-2473, doi: 10.1038/nphys289

[15]

W. L. Shew, H. Yang, S. Yu, R. Roy, and D. Plenz, . Information capacity and transmission are maximized in balanced cortical networks with neuronal avalanches, J. Neurosci.31, 55 (2011).JNRSDS0270-6474

[16]

N. Bertschinger and T. Natschläger, . Real-time computation at the edge of chaos in recurrent neural networks, Neural Comput.16, 1413 (2004).NEUCEB0899-7667

[17]

C. Haldeman and J. M. Beggs, . Critical Branching Captures Activity in Living Neural Networks and Maximizes the Number of Metastable States, Phys. Rev. Lett.94, 058101 (2005).PRLTAO0031-9007

[18]

T. Tanaka, T. Kaneko, and T. Aoyagi, . Recurrent infomax generates cell assemblies, neuronal avalanches, and simple cell-like selectivity, Neural Comput.21, 1038 (2009).NEUCEB0899-7667

[19]

M. Rubinov, O. Sporns, J.-P. Thivierge, and M. Breakspear, . Neurobiologically realistic determinants of self-organized criticality in networks of spiking neurons, PLoS Comput. Biol.7, e1002038 (2011).1553-7358, doi: 10.1371/journal.pcbi.1002038

[20]

J. P. Sethna, K. A. Dahmen, and C. R. Myers, . Crackling noise, Nature (London)410, 242 (2001).NATUAS0028-0836

[21]

J. P. Gleeson and R. Durrett, . Temporal profiles of avalanches on networks, Nat. Commun.8, 1227 (2017).2041-1723, doi: 10.1038/s41467-017-01212-0

[22]

E. J. Friedman and A. S. Landsberg, . Hierarchical networks, power laws, and neuronal avalanches, Chaos23, 013135 (2013).CHAOEH1054-1500

[23]

M. E. J. Newman, . Power laws, Pareto distributions, and Zipf's law, Contemp. Phys.46, 323 (2005).CTPHAF0010-7514

[24]

M. P. H. Stumpf and M. A. Porter, . Critical truths about power laws, Science335, 665 (2012).SCIEAS0036-8075

[25]

M. Martinello, J. Hidalgo, A. Maritan, S. di Santo, D. Plenz, and M. A. Muñoz, . Neutral Theory and Scale-Free Neural Dynamics, Phys. Rev. X7, 041071 (2017)., doi: 10.1103/PhysRevX.7.041071

[26]

M. Benayoun, J. D. Cowan, W. van Drongelen, and E. Wallace, . Avalanches in a stochastic model of spiking neurons, PLOS Comput. Biol.6, e1000846 (2010).1553-7358, doi: 10.1371/journal.pcbi.1000846

[27]

J. Touboul and A. Destexhe, . Power-law statistics and universal scaling in the absence of criticality, Phys. Rev. E95, 012413 (2017).2470-0045, doi: 10.1103/PhysRevE.95.012413

[28]

J. Touboul and A. Destexhe, . Can power-law scaling and neuronal avalanches arise from stochastic dynamics?PLoS One5, e8982 (2010).1932-6203, doi: 10.1371/journal.pone.0008982

[29]

V. Priesemann and O. Shriki, . Can a time varying external drive give rise to apparent criticality in neural systems?PLoS Comput. Biol.14, e1006081 (2018).1553-7358, doi: 10.1371/journal.pcbi.1006081

[30]

M. E. J. Newman, Networks: An Introduction (Oxford University Press, Oxford, 2010).

[31]

K. B. Athreya and P. E. Ney, Branching Processes (Springer-Verlag, Berlin, 1972).

[32]

K.-I. Goh, D.-S. Lee, B. Kahng, and D. Kim, . Sandpile on Scale-Free Networks, Phys. Rev. Lett.91, 148701 (2003).PRLTAO0031-9007

[33]

J. P. Gleeson, J. A. Ward, K. P. O'Sullivan, and W. T. Lee, . Competition-Induced Criticality in a Model of Meme Popularity, Phys. Rev. Lett.112, 048701 (2014).PRLTAO0031-9007

[34]

A.-L. Barabási, H. Jeong, Z. Néda, E. Ravasz, A. Schubert, and T. Vicsek, . Evolution of the social network of scientific collaborations, Phys. A (Amsterdam, Neth.)311, 590 (2002).PHYADX0378-4371

[35]

M. Shimono and J. M. Beggs, . Functional clusters, hubs, and communities in the cortical microconnectome, Cereb. Cortex25, 3743 (2014)., doi: 10.1093/cercor/bhu252

[36]

A.-L. Barabási, . Scale-free networks: A decade and beyond, Science325, 412 (2009).SCIEAS0036-8075

[37]

A. D. Broido and A. Clauset, . Scale-free networks are rare, Nat. Commun.10, 1017 (2019)., doi: 10.1038/s41467-019-08746-5

[38]

J. Goldenberg, B. Libai, and E. Muller, . Talk of the network: A complex systems look at the underlying process of word-of-mouth, Marketing Lett.12, 211 (2001).0923-0645, doi: 10.1023/A:1011122126881

[39]

J. Goldenberg, B. Libai, and E. Muller, . Using complex systems analysis to advance marketing theory development: Modeling heterogeneity effects on new product growth through stochastic cellular automata, Acad. Marketing Sci. Rev.9, 1 (2001).

[40]

D. Kempe, J. Kleinberg, and É. Tardos, Maximizing the spread of influence through a social network, in Proceedings of the Ninth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (ACM, New York, 2003), pp. 137–146.

[41]

A. N. Burkitt, . A review of the integrate-and-fire neuron model: I. Homogeneous synaptic input, Biol. Cybern.95, 1 (2006).BICYAF0340-1200

[42]

A. N. Burkitt, . A review of the integrate-and-fire neuron model: II. Inhomogeneous synaptic input and network properties, Biol. Cybern.95, 97 (2006).BICYAF0340-1200

[43]

A. Levina, J. M. Herrmann, and T. Geisel, . Dynamical synapses causing self-organized criticality in neural networks, Nat. Phys.3, 857 (2007).1745-2473, doi: 10.1038/nphys758

[44]

A. Faqeeh, Percolation and its relations to other processes in networks, Ph.D. thesis, University of Limerick, 2016.

[45]

E. A. Bender and E. Canfield, . The asymptotic number of labeled graphs with given degree sequences, J. Comb. Theory, Ser. A24, 296 (1978).JCBTA70097-3165

[46]

S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D.-U. Hwang, . Complex networks: Structure and dynamics, Phys. Rep.424, 175 (2006).PRPLCM0370-1573

[47]

M. Molloy and B. Reed, . A critical point for random graphs with a given degree sequence, Random Struct, Algor.6, 161 (1995).1042-9832, doi: 10.1002/rsa.3240060204

[48]

See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevE.100.010401 for more details of calculations and/or further relevant results.

[49]

M. E. J. Newman, S. H. Strogatz, and D. J. Watts, . Random graphs with arbitrary degree distributions and their applications, Phys. Rev. E64, 026118 (2001).1063-651X, doi: 10.1103/PhysRevE.64.026118

[50]

M. E. J. Newman, . Spread of epidemic disease on networks, Phys. Rev. E66, 016128 (2002).1063-651X, doi: 10.1103/PhysRevE.66.016128

[51]

R. Cohen, D. ben-Avraham, and S. Havlin, . Percolation critical exponents in scale-free networks, Phys. Rev. E66, 036113 (2002).1063-651X, doi: 10.1103/PhysRevE.66.036113

[52]

N. Schwartz, R. Cohen, D. ben-Avraham, A.-L. Barabási, and S. Havlin, . Percolation in directed scale-free networks, Phys. Rev. E66, 015104(R) (2002).1063-651X, doi: 10.1103/PhysRevE.66.015104

[53]

D. Lee, J. Kim, B. Kahng, and D. Kim, . Scale-free random branching trees in supercritical phase, J. Phys. A: Math. Theor.40, 7139 (2007).1751-8113, doi: 10.1088/1751-8113/40/26/002

[54]

H. S. Wilf, Generating Functionology (A K Peters, Ltd., Wellesley, MA, 2006).

[55]

G. Weiss, Aspects and Applications of the Random Walk, Random Materials and Processes (North-Holland, New York, 2005).

[56]

D. S. Callaway, M. E. J. Newman, S. H. Strogatz, and D. J. Watts, . Network Robustness and Fragility: Percolation on Random Graphs, Phys. Rev. Lett.85, 5468 (2000).PRLTAO0031-9007

[57]

M. E. J. Newman, . Component sizes in networks with arbitrary degree distributions, Phys. Rev. E76, 045101(R) (2007).PLEEE81539-3755

[58]

In fact the property $\stackrel{~}{G}\left(z\right)={G}_{n}\left(1-p+p+p\phantom{\rule{0.16em}{0ex}}z\right)$ (that holds for site percolation too) is used in Refs. [51,52] only to calculate the location of the critical point of site percolation and not to construct the governing equations.

[59]

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, we will later show that in some cases the prefactors can indicate the range of validity of the theoretical predictions.

[60]

P. D. Miller, Applied Asymptotic Analysis, Graduate studies in mathematics (American Mathematical Society, Providence, RI, 2006).

[61]

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).

[62]

G. W. Gross, B. K. Rhoades, H. M. Azzazy, and M.-C. Wu, . The use of neuronal networks on multielectrode arrays as biosensors, Biosens. Bioelectron.10, 553 (1995).BBIOE40956-5663

[63]

J. Müller, M. Ballini, P. Livi, Y. Chen, M. Radivojevic, A. Shadmani, V. Viswam, I. L. Jones, M. Fiscella, R. Diggelmann, A. Stettler, U. Frey, D. J. Bakkum, and A. Hierlemann, . High-resolution CMOS MEA platform to study neurons at subcellular, cellular, and network levels, Lab Chip15, 2767 (2015).1473-0197, doi: 10.1039/C5LC00133A

[64]

T. L. Ribeiro, M. Copelli, F. Caixeta, H. Belchior, D. R. Chialvo, M. A. L. Nicolelis, and S. Ribeiro, . Spike avalanches exhibit universal dynamics across the sleep-wake cycle, PLoS One5, e14129 (2010)., doi: 10.1371/journal.pone.0014129

[65]

E. Tagliazucchi, P. Balenzuela, D. Fraiman, and D. Chialvo, . Criticality in large-scale brain FMRI dynamics unveiled by a novel point process analysis, Front. Physiol.3, 15 (2012).1664-042X, doi: 10.3389/fphys.2012.00015

[66]

A. Faqeeh, S. Melnik, P. Colomer-de-Simón, and J. P. Gleeson, . Emergence of coexisting percolating clusters in networks, Phys. Rev. E93, 062308 (2016).2470-0045, doi: 10.1103/PhysRevE.93.062308

[67]

A. Faqeeh, S. Melnik, and J. P. Gleeson, . Network cloning unfolds the effect of clustering on dynamical processes, Phys. Rev. E91, 052807 (2015).PLEEE81539-3755

[68]

M. A. Muñoz, R. Juhász, C. Castellano, and G. Ódor, . Griffiths Phases on Complex Networks, Phys. Rev. Lett.105, 128701 (2010).PRLTAO0031-9007

[69]

S. Squires, E. Ott, and M. Girvan, . Dynamical Instability in Boolean Networks as a Percolation Problem, Phys. Rev. Lett.109, 085701 (2012).PRLTAO0031-9007

Citing articles via
https://www.researchpad.co/tools/openurl?pubtype=article&doi=10.1103/PhysRevE.100.010401&title=Emergence of power laws in noncritical neuronal systems&author=Ali Faqeeh,Saeed Osat,Filippo Radicchi,James P. Gleeson,&keyword=&subject=Rapid Communications,Biological Physics,