Physical Review. E

American Physical Society

Emergence of power laws in noncritical neuronal systems

Volume:
100,
Issue: 1

DOI 10.1103/PhysRevE.100.010401

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.

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

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

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}/\langle k\rangle $ for UNs or ${q}_{jk}=j{P}_{jk}/\langle j\rangle $ 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 ${\widehat{G}}_{n}\left(z\right)={G}_{n}(1-p+p\phantom{\rule{0.16em}{0ex}}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
$\begin{array}{c}\hfill {H}_{1}\left(z\right)=z{G}_{1}[1-p+p{H}_{1}\left(z\right)],\end{array}$

$\begin{array}{c}\hfill {H}_{0}\left(z\right)=z{G}_{0}[1-p+p{H}_{1}\left(z\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].

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}(1-w)\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}(1-p\varphi )\sim {G}_{1}\left(1\right)-{G}_{1}^{\prime}\left(1\right)p\varphi +{G}_{1}^{\u2033}\left(1\right){\left(p\varphi \right)}^{2}+D{\varphi}^{\stackrel{~}{\lambda}-1}\hfill \\ & \phantom{\rule{1em}{0ex}}+o({\varphi}^{2},{\varphi}^{\stackrel{~}{\lambda}-1})=1-\frac{p}{{p}_{\mathrm{c}}}\varphi +B{\varphi}^{2}+D{\varphi}^{\stackrel{~}{\lambda}-1}\hfill \\ & \phantom{\rule{83.5pt}{0ex}}+o({\varphi}^{2},{\varphi}^{\stackrel{~}{\lambda}-1}),\hfill \end{array}$

$\begin{array}{cc}& {G}_{0}(1-p\varphi )\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}}-{\displaystyle \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}(1-w)& \doteq 1-\varphi \sim (1-w)\left(1-\frac{p}{{p}_{\mathrm{c}}}\varphi +B{\varphi}^{2}+D{\varphi}^{\stackrel{~}{\lambda}-1}\right)\hfill \\ & \Rightarrow \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}(1-w)\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}^{-[1+1/(\stackrel{~}{\lambda}-1)]}$ 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}\langle k\rangle \langle {k}^{3}\rangle}{{\delta}^{2}{\langle {k}^{2}\rangle}^{2}}$ for UNs and ${s}^{*}\approx \frac{2{p}^{2}\langle k\rangle \langle j{k}^{2}\rangle}{{\delta}^{2}{\langle jk\rangle}^{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 \widehat{R}\left(1+\frac{p{p}_{\mathrm{c}}\langle k\rangle}{-\delta}\right){\left(\frac{p{p}_{\mathrm{c}}}{-\delta}\right)}^{{\lambda}_{\mathrm{o}}-1}{s}^{-{\lambda}_{\mathrm{o}}}& \hfill \text{(UDN)}\\ \\ \hfill \widehat{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 $\widehat{R}={a}^{*}(1-\overline{u})$, 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}$.
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$.

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(1-\stackrel{~}{h})=1-p+p\stackrel{~}{h}$. Now, we can write
${H}_{1}\left(z\right)=z\phantom{\rule{0.16em}{0ex}}{G}_{1}(\overline{\eta}+p\epsilon )$

$=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}}{\textstyle \frac{1}{2}}{G}_{1}^{\u2033}\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}}[{a}_{0}+{a}_{1}\epsilon +{a}_{2}{\epsilon}^{2}+\cdots ],$

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}^{\u2033}\left(\overline{\eta}\right){p}^{2}$. Equation (15) gives
$\epsilon \sim \frac{1-{a}_{1}z}{2{a}_{2}z}\pm \frac{1}{2{a}_{2}z}\sqrt{{(1-{a}_{1}z)}^{2}+4{a}_{0}{a}_{2}z(1-z)}.$

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}}\pm \frac{1-{a}_{1}}{2{a}_{2}}\sqrt{1+{s}^{*}(1-z)},$

where ${s}^{*}=\frac{4{a}_{0}{a}_{2}}{{(1-{a}_{1})}^{2}}$. Now,
${H}_{0}\left(z\right)=z\phantom{\rule{0.16em}{0ex}}{G}_{0}(\overline{\eta}+p\epsilon )$

$\sim (1-w)[{G}_{0}\left(\overline{\eta}\right)+{G}_{0}^{\prime}\left(\overline{\eta}\right)p\epsilon +o\left({\epsilon}^{2}\right)]$

$\sim \text{analytical}\phantom{\rule{4.pt}{0ex}}\text{terms}\pm b\sqrt{1+{s}^{*}\left(1-z\right)},$

where $b=\frac{p{G}_{0}^{\prime}\left(\overline{\eta}\right)\phantom{\rule{0.16em}{0ex}}(1-{a}_{1})}{2{a}_{2}}$. Then, according to the asymptotic properties of generating functions [48,54,55],
${\pi}_{s}\sim \pm \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 \pm {\displaystyle \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}^{\u2033}\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}^{\u2033}\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)].
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).

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

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

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

Citing articles via

Tweets

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,

© 2020 Newgen KnowledgeWorks | Privacy & Cookie Policy | Powered by: Nova