Home | About | Journals | Submit | Contact Us | Français |

**|**Scientific Reports**|**PMC5541036

Formats

Article sections

Authors

Related links

Sci Rep. 2017; 7: 7107.

Published online 2017 August 2. doi: 10.1038/s41598-017-07135-6

PMCID: PMC5541036

Afshin Montakhab, Email: ri.ca.uzarihs@bahkatnom.

Received 2016 December 2; Accepted 2017 June 23.

Copyright © The Author(s) 2017

Networks of excitable nodes have recently attracted much attention particularly in regards to neuronal dynamics, where criticality has been argued to be a fundamental property. Refractory behavior, which limits the excitability of neurons is thought to be an important dynamical property. We therefore consider a simple model of excitable nodes which is known to exhibit a transition to instability at a critical point (*λ*=1), and introduce refractory period into its dynamics. We use mean-field analytical calculations as well as numerical simulations to calculate the activity dependent branching ratio that is useful to characterize the behavior of critical systems. We also define avalanches and calculate probability distribution of their size and duration. We find that in the presence of refractory period the dynamics stabilizes while various parameter regimes become accessible. A sub-critical regime with *λ*<1.0, a standard critical behavior with exponents close to critical branching process for *λ*=1, a regime with 1<*λ*<2 that exhibits an interesting scaling behavior, and an oscillating regime with *λ*>2.0. We have therefore shown that refractory behavior leads to a wide range of scaling as well as periodic behavior which are relevant to real neuronal dynamics.

Various physical, biological and chemical systems are composed of interacting excitable agents and thus networks of excitable nodes are widely used to model the behavior of such systems. Examples of such systems include tectonic plates^{1, 2}, Neural networks^{3–7}, models of self-organized criticality (SOC)^{8–11}, and epidemic (contagion) spreading^{12–16}. Oftentimes excitable nodes are modeled with threshold dynamics as in the case of sandpile models of SOC which are thought to underlie the wide range of scale-invariant behavior seen in Nature.

Criticality in cortical dynamics is by now a widely studied and well established field^{17}. Neuronal avalanches have been reported as collective scale-invariant behavior of neurons in the cortical layers of mammalian brain^{3, 17–25}. Probability distribution functions of duration and size of the neuronal avalanches are power law functions observed in a wide range of space and time which are thought to be hallmarks of critical systems. In addition to avalanche statistics, branching ratio has also been used to characterize various time-series in order to establish critical dynamics of the brain^{3}. Criticality of the brain is therefore the subject of a myriad of theoretical as well as experimental studies^{26–34}. Critical dynamics is believed to underlie many functional advantages in a healthy brain, including learning^{35}, optimal dynamic range^{5, 36–39}, efficient information processing^{32}, as well as optimal transmission and storage of information^{4}.

Many excitable agents often display refractory behavior. This behavior which is characterized by a time scale (i.e. refractory period) is a time during which the excited agent cannot be re-excited. The presence of such refractory period can clearly affect the collective dynamics of excitable nodes. Neuronal systems are perfect examples of networks of excitable nodes with refractory period^{40}.

A particularly useful model of excitable agents considers a random network of such nodes with quenched excitatory connection weights and probabilistic dynamics^{41}. It is well-known that such a model exhibits a phase transition between stable and unstable regimes as the average weight of connections is increased^{37, 38, 42}. More recently, it has been shown that the addition of inhibitory connections leads to a ceaseless dynamics which exhibits critical behavior by fine-tuning the system to the transition point associated with the mentioned instability^{43}. This fine-tuning essentially leads to an intricate balance between excitatory and inhibitory connections which may not be *a priori* available and thus a non-generic behavior. Here, we propose to study the original model (without inhibitory connections) in the presence of refractory period. Interestingly, we find that refractory period leads to stable dynamics in the entire range of parameter. However, and perhaps more interestingly, we are able to identify a critical point with robust finite-size scaling behavior, a wide critical-like regime with interesting scaling behavior, as well as a regime with periodic behavior. Therefore, we show that refractory period in addition to stabilizing the dynamics leads to a wide range of parameters which show scaling or periodic behavior which are hallmarks of real neuronal dynamics.

This article is organized as follows: The following section discusses the model. Next, we show our analytical and numerical results, respectively. We close by providing concluding remarks.

The model consists of *N* excitable nodes on a random directed graph where every two nodes are connected with probability *q*. The average out-degree *and* in-degree of the network is equal to *k*=*qN*. Weights of the connections (*w*
_{ij}), that form the adjacency matrix of the network, are randomly chosen real numbers in the range of [0, 2*σ*] with the average connection weight of *σ*. If node *i* is not connected to node *j*, their connection weight is set to zero (*w*
_{ij}=0). Every node can be in one of active or quiescent states. Activity of the network is shown by the spatio-temporal binary variable, *A*
_{i}(*t*), i.e. if the node *i* is active at time *t* then *A*
_{i}(*t*)=1 and when the node is quiescent *A*
_{i}(*t*)=0. The probability of a node to be activated at time *t*+1 is equal to

$$p({A}_{i}(t+\mathrm{1)}=\mathrm{1)}={\delta}_{\mathrm{0,}{A}_{i}(t)}f\left(\sum _{j=1}^{N}\phantom{\rule{.1em}{0ex}}{w}_{ij}{A}_{j}(t)\right)$$

1

where ${\delta}_{\mathrm{0,}{A}_{i}(t)}$ is the Kronecker delta, which is equal to zero (one) if *A*
_{i}(*t*)=1 (0), implying refractory period of one time step, i.e. the activation probability of a node that is active at time *t* will be equal to zero at *t*+1. *f* is a transfer function that transfers the total input of a node $\left({\sum}_{j=1}^{N}\phantom{\rule{.1em}{0ex}}{w}_{ij}{A}_{j}(t)\right)$ at time *t* to the activation probability of that node at time *t*+1, and is chosen to be

$$f(y)=\{\begin{array}{ll}y,\hfill & 0\le y\le 1\hfill \\ 1,\hfill & y>1.\hfill \end{array}$$

2

In this work, our focus is on the aggregate activity of the system which is defined as the number of active nodes at each time step and is equal to ${x}_{t}={\sum}_{i=1}^{N}\phantom{\rule{.1em}{0ex}}{A}_{i}(t)$. We use *x*
_{t} as a dynamical variable to analyze stability as well as statistical properties of the system.

Before presenting our results, it is instructive to remark on properties of the model without a refractory period. It is well known that the behavior of this model, without a refractory period, is governed by the largest eigenvalue of the adjacency matrix^{37–39, 42} which is equal to *λ*=*σ**k* for the random graph that we use^{42, 44}. It has been shown that this system exhibits stability and instability in activity for *λ*≤1 and *λ*>1, respectively. In the case of stability, the system requires external drive to be activated and *x*=0 is the stable attractor of the dynamics. When unstable (*λ*>1) the activity of the system increases and saturates at *x*=*N*. The critical point of the system is *λ*=1 where the system undergoes a transition from stability to instability. Poised at the critical point, the system exhibits scaling behavior for statistical properties of cascades of activity (avalanches) that start by an external perturbation. Our main goal in this paper is to scrutinize the behavior of the system when refractory period is introduced into the dynamics, i.e. delta function in Eq. (^{1}). We are particularly interested in stability of dynamics and/or whether generic scaling behavior similar to cortical samples could arise in the model.

Dynamics of the system with refractory period is governed by Eq. (^{1}). It is clear that the dynamics is strongly affected by the interaction weights (*w*
_{ij}). We can roughly explain the behavior of the system by considering different extremes of *λ* which is proportional to the average weight of connections *σ*. For small values of *λ* ≪ 1 the interaction weights are small and any activity that starts by an initial perturbation is expected to die out very fast leading to a stability of the fixed point *x*=0. In the other extreme, for large values of *λ* ≫ 1, due to large values of interaction weights the probability of firing is expected to be equal to one for any node that receives input and is also quiescent at the time (see Eq. (^{1})). Therefore, at any given time, after a transient period, there are two sets of nodes: *x* that are active, and *N*−*x* that are in refractory period and will be active in next time step. Clearly, the system exhibits oscillations in this limit where all active nodes become quiescent and vice versa. We generally expect that there exists an important range of *λ* over which the system changes its behavior from a ceasing stable phase to a cease-less periodic phase. It is possible that the transition passes though a critical point or region.

In order to analyze the behavior of the model we use the aggregate activity of the system *x*
_{t} as a function of time, and calculate the *activity dependent* branching ratio *b*(*M*)^{45} which is the expectation value of *x*
_{t+1}/*M* when there are *M* active nodes at time *t*,

$$b(M)=E\left(\frac{{x}_{t+1}}{M}|{x}_{t}=M\right).$$

3

It is clear from definition of *b* that for a given value of *x*
_{t}=*M*, if *b*(*M*)>1 (*b*(*M*)<1) an average increase (decrease) in activity is expected in the next time step. Therefore, *b*(*M*) can provide important statistical information about the behavior of the system for the entire range of possible values of *x*
_{t}.

Galton-Watson theory of branching process holds that a branching ratio less, equal or larger than one respectively bespeaks sub-critical, critical and super-critical phases in a system^{46}. But, the activity dependent branching ratio, due to its variability with respect to activity *x*, is different from the branching ratio defined in the Galton-Watson process and therefore provides much more information about the dynamics of a system. We thus consider criticality with regard to the activity dependent branching ratio. We define a system as critical if there exists a range *R* (*M**R*) which is accessible by the long term dynamics of the system and exhibits two characteristics: (i) the value of the activity dependent branching ratio must be equal to one over *R* in the thermodynamic limit, i.e. lim_{N→∞} *b*(*M* ∈ *R*) = 1, and (ii) *R* must go to infinity as *N*→∞. Condition (i) has to do with critical systems being (on the average) unpredictable. Condition (ii) has to do with lack of characteristic scale for a critical systems in the thermodynamic limit. We note that *b*(*M*) has been used to ascertain criticality in a wide range of systems including sandpile models of SOC or solar flares^{45} as well as neural networks^{33}.

We begin by employing a mean-field approach in order to calculate the activity dependent branching ratio, analytically. If there are *M* active nodes at a time *t*, then at time *t*+1 every node will receive an input with the weight of $\frac{M\langle k\rangle}{N}$, and if a node is quiescent at time *t* will be activated with probability of $f\left(\frac{M\langle k\rangle \sigma}{N}\right)$. The largest eigenvalue of the connection matrix is *λ*=*σ**k* and the activation probability can thus be written as $f\left(\frac{M\lambda}{N}\right)$. Since $\frac{M\lambda}{N}\ge 0$, two situations are possible: (a) $0\le \frac{M\lambda}{N}\le 1$: in this case $f\left(\frac{M\lambda}{N}\right)=\frac{M\lambda}{N}$, and the probability of having exactly *x*
_{t+1}=*z* active nodes at *t*+1 when there are *x*
_{t}=*M* active nodes at time *t* can be approximated as a binomial probability function, i.e. due to the refractory period of one time step, there are *N*−*M* nodes that can be activated with probability $\frac{M\lambda}{N}$. (b) $\frac{M\lambda}{N}>1$: in this case $f\left(\frac{M\lambda}{N}\right)=1$ and every quiescent node that receives input in a time step will be activated in the next time step. Therefore, the conditional probability is obtained as:

$$P({x}_{t+1}=z|{x}_{t}=M)=\{\begin{array}{ll}\left(\begin{array}{c}N-M\\ z\end{array}\right)\phantom{\rule{.25em}{0ex}}{\left(\frac{M\lambda}{N}\right)}^{z}\phantom{\rule{.25em}{0ex}}{\left(1-\frac{M\lambda}{N}\right)}^{N-M-z}\hfill & \frac{M\lambda}{N}<1\hfill \\ {\delta}_{z,N-M}\hfill & \frac{M\lambda}{N}\ge 1\hfill \end{array}$$

4

where the Kronecker delta indicates that all quiescent nodes will on the average be activated when *Mλ*/*N*>1. We can therefore calculate the activity dependent branching ratio as:

$$b(M)=\frac{E({x}_{t+1}|{x}_{t}=M)}{M}=\{\begin{array}{ll}\lambda -\frac{\lambda}{N}M\hfill & M<\frac{N}{\lambda}\hfill \\ \frac{N}{M}-1\hfill & M\ge \frac{N}{\lambda}\hfill \end{array}$$

5

Having calculated the function *b*(*M*), we can present a mean-field analysis of the behavior of the system. As is clear from Eq. (^{5}) the behavior crucially depends on the value of *λ*. Note that *b*(*M*) is a piecewise differentiable decreasing function of *M*, which has a linear as well as a nonlinear regime and goes to zero at *M*=*N*, see Fig. 1. For *λ*<1 (see Fig. 1(a)), we have $M<\frac{N}{\lambda}$ and we are in the linear branch of $b(M)=\lambda -\frac{\lambda}{N}M$ which is always less than one. This indicates that the activity of the system will on average decrease until reaching the fixed point of *x*
_{t}=*x**=0 for all initial perturbations. This is the stable sub-critical regime for which *b*(*M*) < 1, ∀ *M*. Note that, when *λ*→1, *b*(*M*)|_{M=x⁎=0} → 1 and the stable fixed point is expected to exhibit critical behavior. Also note that in the *N*→∞ limit, *b*(*M*) = *λ*, ∀ *M* for *λ*≤1 which indicates subcritical behavior for *λ*<1 *and* critical behavior for *λ*=1.

Activity dependent branching ratio for different parameter regimes of (**a**) *λ*<1, (**b**) 1<*λ*<2, (**c**) *λ*=2 and (**d**) *λ*>2. **...**

For 1<*λ*<2 (see Fig. 1(b)), the line *b*(*M*)=1 will intersect *b*(*M*) in the linear regime at ${M}_{c}=N\left(1-\frac{1}{\lambda}\right)$. Interestingly, due to the negative slope of *b*(*M*), *x**=*M*
_{c} will be the attractor of the dynamics as *M*<*M*
_{c} (*M*>*M*
_{c}) *b*(*M*) will be larger (smaller) than one which would indicate increased (decreased) average activity until *x*=*M*
_{c} is reached where *b*(*M*)=1. This indicates a stable, ceaseless (*x*^{⁎} ≠ 0) dynamics which would exhibit critical behavior as *b*(*M*=*M*
_{c})=1. Note that for *λ*=2 (see Fig. 1(c)), the critical attractor *M*
_{c} will coincide with nonlinear regime of *b*(*M*) at ${M}_{c}=\frac{N}{2}$.

As *λ* is further increased, *λ*>2 (see Fig. 1(d)) the line *b*(*M*)=1 will intersect *b*(*M*) in the nonlinear regime and the simple analysis presented above will no longer hold. However, one can simply note that for $M>\frac{N}{\lambda}$, *E*(*x*_{t+1}|*x*_{t} = *M*) = *N* − *M* (see Eq. (^{5})) and *E*(*x*_{t+1}|*x*_{t} = *N* − *M*) = *M* which would indicate a period-2 oscillating behavior. For the case when *x*
_{t}<*N*/2 (i.e. initial conditions in the linear regime), time evolution of the system will increase *x*
_{t} as *b*(*M*)>1 in this regime until we reach $M=\frac{N}{\lambda}$ after which the same periodic behavior would occur between $M=\frac{N}{\lambda}$ and $M=N-\frac{N}{\lambda}$. The fact that periodic behavior arises as a result of refractory period has previously been observed as in models of epidemic spreading^{47, 48}. However, the case of refractory period larger than one presents an interesting case study, the details of which is presented in the Appendix.

Figure 2 summarizes the mean-field behavior of the system. Three parameter regimes are presented. For 0<*λ*<1 the activity dependent branching ratio at the stable attractor of the dynamics is equal to *b*(*x**)=*λ*<1 and the system is sub critical. For 1.0≤*λ*≤2.0 the activity dependent branching ratio is *b*(*x**)=1.0 at the stable attractor of the dynamics and criticality is expected. For *λ*>2 we have two values for the branching function and the dynamics jumps back and forth between these two values. *λ*
_{c}=1 and *λ*
_{p}=2 are the bifurcation points where the system changes behavior. while the value of *λ*
_{c} is independent of refractory period, the value of *λ*
_{p} (>1) depends on the choice of refractory period, see Appendix.

The above results portray the general behavior of the system in mean-field approximation where fluctuations have been ignored. However, as is well-known fluctuations are very important and tend to dominate system behavior in the critical regime. We therefore propose to study our system by extensive numerical simulations for *N*=1×10^{4} up to 8×10^{4} and *q*=0.01 in the range of 0.9≤*λ*<2 with particular focus on the critical behavior of the system. All activities are initiated by choosing a random site *i* and setting *A*
_{i}(*t*=0)=1 and following the ensuing dynamics according to Eqs (^{1} and ^{2}). *x*
_{t} is recorded for long times from which we can easily calculate *b*(*M*) numerically.

Figure 3 shows our results for *b*(*M*) for various system sizes and *λ*=1.2 as well as *λ*=1.0, where criticality is expected. In both cases, our results clearly show a linear behavior in accordance with Eq. (^{5}), where $b(M)=\lambda -\frac{\lambda}{N}M$ provides a prefect fit to the data. Note that as *N* is increased the range of system’s activity increases as the slope ($\frac{\lambda}{N}$) goes to zero. Therefore, one can expect that in the large system size limit *b*(*M*)→1 for all accessible *M* indicating a critical behavior. We have also checked various other values of *λ* and have observed similar behavior to that of *λ*=1.2 (above) for the range of 1<*λ*<2 (not shown).

In order to better understand the behavior of fluctuations about the attractors *x*
^{*}, we have plotted the probability distribution function of system’s activity as *p*(*x*) in Fig. 4 for *λ*=1.2, 1.0, and 0.9. For *λ*=1.2 (Fig. 4(a)) we observe a Gaussian behavior which peaks exactly at $x={M}_{c}=N\left(1-\frac{1}{\lambda}\right)$. It is interesting to note that our results indicate that the width of the Gaussian increases with the system size as *ξ* ∼ *N*^{0.5} (see inset of Fig. 4(a)), in accordance with the central limit theorem. For *λ*=1.0 (Fig. 4(b)), we observe a distinctly different behavior as *p*(*x*) displays a power-law behavior with system size dependent cutoff. It is, however, important to note that *p*(*x*) is maximized at *x*=0 as indicated by our mean-field analysis. In Fig. 4(c), we plot the same results for *λ*=0.9 for various system sizes on a log-linear plot, all of which coincide on the same curve. We therefore conclude that *p*(*x*) displays an exponential behavior with a scale which is system size *independent*. Again, the attractor *x*=0 appears as the most probable state, however, the size independent scale in (c) as opposed to size dependent scale in part (b) (and even (a)) is the distinction between sub-critical and critical systems.

Probability distribution function of aggregate activity (*p*(*x*)) as a function of *x* for different values of *N*. Panel (a) is a linear plot showing a Gaussian function for *λ*=1.2 with a maximum at *x*=*x*
_{c}= **...**

As we have defined a critical system based on the behavior of *b*(*M*) (see above), we have shown that our system exhibits critical behavior in the range 1≤*λ*<2. However, distinctly different behavior is observed for *P*(*x*) as *M*
_{c}=0 (Fig. 4(b)) changes to *M*_{c} ≠ 0 (Fig. 4(a)). To better understand the critical behavior of the system we now turn our attention to avalanches. In the case of *λ*=1.0 for which the stable attractor of the dynamics is *x**=0, the avalanches are well defined as the activity between two stabilities initiated by an external perturbation. For other values of 1<*λ*<2 over which the system exhibits self-sustaining behavior we define the avalanches (see Fig. 5), as the continuous aggregate activity of nodes above a threshold value *x*
_{th}. The number of time steps of an excursion above *x*
_{th} is defined as duration (D) and the summation ∑_{D} *x*_{t} − *x*_{th} as the size (*S*) of an avalanche.

An avalanche is defined as excursion of the aggregate activity of nodes above a threshold value *x*
_{th}, *D* is the defined as the duration and the integral between *x*
_{t} and *x*
_{th} (the colored area) as the size (*S*) of the avalanche.

Probability distribution function of size and duration of avalanches (*P*(*y*), *y* ∈ {*S*, *D*}) are calculated for systems with different *N* over the parameter regime 1≤*λ*≤2. In the case of criticality, these probability distribution functions are expected to exhibit power-law behavior with a cutoff which is an increasing function of *N*. The usual scaling ansatz for such a behavior is *P*(*y*) ∼ *y*^{−τy}*g*_{y}(*y*/*N*^{βy}), where *g*
_{y} is a universal cutoff function that is identical for different system sizes. *τ*
_{y} is the critical exponent and *β*
_{y} is referred to as the finite-size scaling exponent. When criticality holds, if we rescale *y* → *y*/*N*^{βy} and *P*(*y*) → *y*^{τy}*P*(*y*) then the plots of rescaled variables must collapse into one universal curve for the correct values of *τ*
_{y} and *β*
_{y}
^{11}.

Probability distribution functions of *S* are plotted in the main panel of Fig. 6 for systems with *λ*=1.2 and different sizes. It is clearly seen that the plots exhibit a power-law region and a cutoff that increases by the system size. Inset panel of Fig. 6 shows collapse of the rescaled data of the main panel with exponents *τ*
_{S}=1.00±0.01 and *β*
_{S}=0.50±0.01. It must be noted that our numerical analysis show that different choices of *x*
_{th} does not change the values of *τ*
_{S} and *β*
_{S}, and we choose *x*
_{th}=*x*
_{c} where we have better statistics for avalanche distribution functions. Our numerical analysis for other values of 1<*λ*<2 indicate that the same values of *τ*
_{S}=1.00±0.01 and *β*
_{S}=0.50±0.01 are also obtained.

Main panel: probability distribution functions of avalanche sizes for systems with *λ*=1.2 and different values of *N*. Inset panel: plots of rescaled data, collapsed into one universal curve with *τ*
_{S}=1.0 **...**

However, the critical behavior obtained for *λ*=1.0 is somewhat different. As shown in Fig. 7(a), we obtain *τ*
_{S}=1.46±0.02 and *β*
_{S}=1.00±0.03 which are indication of a different universality class for *λ*=1.0, where the critical exponent is close to the critical branching process, i.e. *τ*
_{S}=3/2. More importantly, however, our study of finite-size scaling of avalanche sizes, in addition to *b*(*M*)→1 presented earlier, provide firm evidence for critical behavior of the model in the range of 1≤*λ*<2. This could potentially provide an explanation to a wide range of criticality observed in neuronal systems, without any apparent tuning of parameters.

Main panels: probability distribution functions of (**a**) avalanche sizes and (**b**) avalanche durations for systems with *λ*=1.0 and different values of *N*. Inset panel: plots of rescaled data, collapsed into one universal curve with **...**

We have also calculated probability distribution functions of avalanche durations (*P*(*D*)). In Fig. 7(b), we have shown reslults for the *λ*=1.0 case, where we obtain *τ*
_{D}=1.86±0.02 and *β*
_{D}=0.50±0.01 again close to the exponent *τ*
_{D}=2.0 of the critical branching process. However, the model exhibits an unusual scaling behavior for *D* in 1<*λ*<2 range. For system sizes we have been able to study (i.e up to *N*=8×10^{4}), probability distribution functions of avalanche durations *do* exhibit a power-law behavior. However, no appreciable increase is observed in the cutoff for the present system sizes. This behavior could possibly happen if the finite-size exponent *β*
_{D} is so small that the cutoff function does not change considerably for the system sizes we have considered here.

To shed light on such a behavior, we consider a scaling anzats that relates the size and duration of avalanches as:

6

in which *E*(*S*|*D*) is the expectation value of *S* when *D* is given. As seen in Fig. 8, *E*(*S*|*D*) is a linear function of *D* for a given system size. Careful regression analysis shows that *α*=0.50±0.01 and *γ*=1.00±0.01. On the other hand, from the above numerical analysis, we know that the maximum value of avalanche sizes scales as *S*_{max} ∼ *N*^{βS}. Due to the linear relation between *E*(*S*|*D*) and *D* we can write *E*(*S*|*D*
_{max})=*S*
_{max}. Using Eq. (^{6}) we write *S*_{max} ∼ *N*^{α}*D*_{max} which leads to *D*_{max} ∼ *N*^{βS−α} and therefore gives *β*
_{D}=*β*
_{S}−*α*=0.00±0.02. The fact that *β*_{D} ≈ 0 indicates why we do not observe finite size scaling for avalanche durations despite the fact that we observe a power-law behavior for *P*(*D*) in a limited range of data. This is an interesting case whose full understanding requires further investigation with much larger system sizes than studied here.

An interesting property of critical dynamical systems, such as self-organized critical models, is that short-term evolution of perturbations is a power-law function of time and the system exhibits power-law sensitivity to initial conditions^{49}. In order to present another evidence for criticality of the system in the parameter regime 1≤*λ*<*λ*
_{p}, we test this behavior in our model. In order to do that we consider a system in this parameter regime and run the simulation until the system reaches its critical” state. At this point we pause the simulation for a moment. We have the activity vector **A**(*t*) = {*A*_{1}(*t*), *A*_{2}(*t*), *A*_{3}(*t*), *…*, *A*_{N}(*t*)} in an *N* dimensional space and make a copy of the system **A**′. Then, we introduce a small perturbation to **A**′ by changing a few randomly chosen elements (${\mathbf{A}}_{i}^{\prime}$) from one to zero or vice versa. The difference between two systems, which is a distance in the N dimensional space, is defined as the Hamming distance (*H*) between **A** and **A**′ which is

$$H(t)=\sum _{1}^{N}\phantom{\rule{.1em}{0ex}}{({A}_{i}(t)-{A}_{i}^{\prime}(t))}^{2}$$

7

Now, we can study the short-term evolution of *H*(*t*). It must be noted that we used the same random seed for simulating both systems. In order to have firm results we need to do an ensemble averaging, therefore, we have done the above process for 2000 realizations. Time evolution of the ensemble-averaged Hamming distance is plotted in Fig. 9 for systems with refractory period of one time step, *N*=40000 and different values of *λ*=1.0, 1.2, 1.4, 1.6, 1.8. It is observed that the system exhibits power-law sensitivity to initial conditions

8

This provides yet another evidence for criticality of the system in the parameter range of 1≤*λ*<*λ*
_{p}. Note that we have also included *λ*
_{c}=1.0 where criticality is well established. The exponent *δ* is an increasing function of *λ* indicating more chaotic behavior.

Motivated by the fact that scaling behavior is observed in systems composed of interacting excitable nodes, and that excitable nodes can exhibit refractory period, as in neuronal systems, we have studied a simple model of excitable nodes on a random directed graph in the presence of refractory period. The behavior of the model without refractory period has been well-studied previously.

We find that in the presence of refractory period the behavior of system undergoes dramatic changes and different dynamical regimes become accessible for the system. A sub-critical regime for *λ*<1.0. A standard critical behavior with scaling and power-law behavior of avalanche sizes and durations, similar to the critical branching process, for *λ*=1. An important regime with stable self-sustaining dynamics with interesting scaling behavior for 1<*λ*<*λ*
_{p} where activity dependent branching ratio goes to one in the thermodynamic limit and the system exhibits power-law statistics in avalanche distribution functions as well as power-law sensitivity to initial perturbations. The critical exponents associated with this new regime are distinct from those obtained from the standard critical point at *λ*=1. Finally an oscillating regime with *λ*>*λ*
_{p} where the dynamics oscillates between various well-defined states is observed. While our results highlights the importance of refractory period in networks of excitable nodes, it can also provide understanding for similar behavior observed in real neuronal systems. Branching ratios equal to one, as well as power law statistics in neuronal avalanches have been observed in a wide variety of real neuronal systems, with no apparent tuning of a parameter. The fact that our model exhibits similar behavior for a wide range of parameter (i.e. 1≤*λ*<*λ*
_{p}) is the main result of our study. The more general case of longer refractory periods does not change our general conclusions and the details of such generalization appear in the Appendix.

In a recent work Larremore *et al*. have shown that inhibition causes ceaseless dynamics at or near the critical point ($\lambda \underset{\u02dc}{<}1$) in a similar model^{43}. We, on the other hand, have also observed ceaseless stable dynamics but in a different parameter regime of *λ*>1, when refractory period is included without inhibition. Both inhibition and refractory period are thought to decrease the level of activity in a system. However, they seem to lead to stable ceaseless dynamics in networks of excitable nodes. It would be interesting to study the effect of both these important properties simultaneously to see if larger parameter regime with stable dynamics and scaling (critical) behavior can be obtained.

It is generally believed that networks with small-world effect would exhibit critical behavior with exponents associated with critical branching process which correspond to mean-field exponents *τ*
_{S}=3/2 and *τ*
_{D}=2. Here, we have observed non-mean-field behavior with exponents which are significantly smaller than the critical branching process for a wide range of parameter in a small-world network of excitable nodes with refractory period. Of course, we have also obtained similar mean-field exponents for a particular value of *λ*=1.

Although neuronal dynamics have been the main motivation of our work, other important dynamical processes such as epidemic spreading could also have relevance to our work as refractory period is thought to be important in such spreading processes as well.

Computer code simulations were developed in FORTRAN 90. In order to get good statistics for probability distribution functions of *x* as well as duration and size of avalanches, for each system size (*N*=10^{4}, 2×10^{4}, 4×10^{4}, 8×10^{4}) we performed simulations for 20 different realization of networks and 10^{6} time steps for each network. Random networks were made by having every two node connected with a probability of *q*=0.01.

Support of Iranian National Elites Foundation, Shiraz University Research Council and Institute for Advanced Studies in Basic sciences (IASBS) is acknowledged. S.A.M. also acknowledges Y. Souboti without whose support this work could not be completed.

Author Contributions

S.A.M. and A.M. conceived the project. S.A.M. carried out the numerical simulations as well as analytical calculations. S.A.M., A.M. and A.V. analyzed the results. S.A.M. and A.M. wrote the manuscript. All authors reviewed and approved the manuscript.

The authors declare that they have no competing interests.

S. Amin Moosavi and Afshin Montakhab contributed equally to this work.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

1. Burridge R, Knopoff L. Model and theoretical seismicity. Bull. Seismol. Soc. Am. 1967;57:341–371.

2. Olami Z, Feder HJS, Christensen K. Self-organized criticality in a continuous, nonconservative cellular automaton modeling earthquakes. Phys. Rev. Lett. 1992;68:1244. doi: 10.1103/PhysRevLett.68.1244. [PubMed] [Cross Ref]

3. Beggs JM, Plenz D. Neuronal avalanches in neocortical circuits. J. Neurosci. 2003;23:11167–11177. [PubMed]

4. Shew WL, Yang H, Yu S, Roy R, Plenz D. Information capacity and transmission are maximized in balanced cortical networks with neuronal avalanches. J. Neurosci. 2011;31:55–63. doi: 10.1523/JNEUROSCI.4637-10.2011. [PMC free article] [PubMed] [Cross Ref]

5. Shew WL, Yang H, Petermann T, Roy R, Plenz D. Neuronal avalanches imply maximum dynamic range in cortical networks at criticality. J. Neurosci. 2009;29:15595–15600. doi: 10.1523/JNEUROSCI.3864-09.2009. [PMC free article] [PubMed] [Cross Ref]

6. Haldeman C, Beggs JM. Critical branching captures activity in living neural networks and maximizes the number of metastable states. Phys. rev. lett. 2005;94:058101. doi: 10.1103/PhysRevLett.94.058101. [PubMed] [Cross Ref]

7. Williams-García RV, Moore M, Beggs JM, Ortiz G. Quasicritical brain dynamics on a nonequilibrium widom line. Phy. Rev. E. 2014;90:062714. doi: 10.1103/PhysRevE.90.062714. [PubMed] [Cross Ref]

8. Bak P, Tang C, Wiesenfeld K. Self-organized criticality: An explanation of the 1/*f* noise. Phys. Rev. Lett. 1987;59:381. doi: 10.1103/PhysRevLett.59.381. [PubMed] [Cross Ref]

9. Bak P, Tang C, Wiesenfeld K. Self-organized criticality. Phys. Rev. A. 1988;38:364. doi: 10.1103/PhysRevA.38.364. [PubMed] [Cross Ref]

10. Bak, P. *How nature works*: *the science of self*-*organized criticality* (Springer-Verlag New York, 1996).

11. Pruessner, G. *Self*-*organized criticality*: *theory*, *models and characterisation* (Cambridge University Press New York, 2012).

12. Son SW, Bizhani G, Christensen C, Grassberger P, Paczuski M. Percolation theory on interdependent networks based on epidemic spreading. Europhys. Lett. 2012;97:16006. doi: 10.1209/0295-5075/97/16006. [Cross Ref]

13. Karrer B, Newman MEJ. Competing epidemics on complex networks. Phys. Rev. E. 2011;84:036106. doi: 10.1103/PhysRevE.84.036106. [PubMed] [Cross Ref]

14. Van Mieghem P. Epidemic phase transition of the SIS type in networks. Erophys. Lett. 2012;97:48004. doi: 10.1209/0295-5075/97/48004. [Cross Ref]

15. Montakhab A, Manshour P. Low prevalence, quasi-stationarity and power-law behavior in a model of contagion spreading. Erophys. Lett. 2012;99:58002. doi: 10.1209/0295-5075/99/58002. [Cross Ref]

16. Montakhab A, Manshour P. Contagion spreading on complex networks with local deterministic dynamics. Commun. Nonlinear Sci. Numer. Simul. 2014;19:2414–2422. doi: 10.1016/j.cnsns.2013.12.015. [Cross Ref]

17. Plenz, D. ed. *Criticality in Neural Systems* (New York, NY: Wiley-VCH, 2014).

18. Plenz D, Thiagarajan TC. The organizing principles of neuronal avalanches: cell assemblies in the cortex? Trends in Neurosci. 2007;30:101–110. doi: 10.1016/j.tins.2007.01.005. [PubMed] [Cross Ref]

19. Beggs JM, Plenz D. Neuronal avalanches are diverse and precise activity patterns that are stable for many hours in cortical slice cultures. J. Neurosci. 2004;24:5216–5229. doi: 10.1523/JNEUROSCI.0540-04.2004. [PubMed] [Cross Ref]

20. Friedman N, et al. Universal critical dynamics in high resolution neuronal avalanche data. Phys. Rev. Lett. 2012;108:208102. doi: 10.1103/PhysRevLett.108.208102. [PubMed] [Cross Ref]

21. Petermann T, et al. Spontaneous cortical activity in awake monkeys composed of neuronal avalanches. Proc. Natl. Acad. Sci. USA. 2009;106:15921–15926. doi: 10.1073/pnas.0904089106. [PubMed] [Cross Ref]

22. Tagliazucchi E, Balenzuela P, Fraiman D, Chialvo DR. Criticality in large-scale brain fMRI dynamics unveiled by a novel point process analysis. Front. Physiol. 2012;3:15. doi: 10.3389/fphys.2012.00015. [PMC free article] [PubMed] [Cross Ref]

23. Shriki O, et al. Neuronal avalanches in the resting MEG of the human brain. J. Neurosci. 2013;33:7079–7090. doi: 10.1523/JNEUROSCI.4286-12.2013. [PMC free article] [PubMed] [Cross Ref]

24. Haimovici A, Tagliazucchi E, Balenzuela P, Chialvo DR. Brain organization into resting state networks emerges at criticality on a model of the human connectome. Phys. Rev. Lett. 2013;110:178101. doi: 10.1103/PhysRevLett.110.178101. [PubMed] [Cross Ref]

25. Chialvo D. Emergent complex neural dynamics. Nature Physics. 2010;6:744–750. doi: 10.1038/nphys1803. [Cross Ref]

26. Beggs JM, Timme N. Being critical of criticality in the brain. Front. in physiol. 2012;3:163. doi: 10.3389/fphys.2012.00163. [PMC free article] [PubMed] [Cross Ref]

27. Chialvo D. Critical brain networks. Physica A. 2004;340:756–765. doi: 10.1016/j.physa.2004.05.064. [Cross Ref]

28. Bornholdt S, Röhl T. Self-organized critical neural networks. Phys. Rev. E. 2003;67:066118. doi: 10.1103/PhysRevE.67.066118. [PubMed] [Cross Ref]

29. Levina A, Herrmann JM, Geisel T. Dynamical synapses causing self-organized criticality in neural networks. Nature Physics. 2007;3:857–860. doi: 10.1038/nphys758. [Cross Ref]

30. Millman D, Mihalas S, Kirkwood A, Niebur E. Self-organized criticality occurs in non-conservative neuronal networks during up’ states. Nature physics. 2010;6:801–805. doi: 10.1038/nphys1757. [PMC free article] [PubMed] [Cross Ref]

31. de Arcangelis L, Perrone-Capano C, Herrmann HJ. Self-Organized Criticality Model for Brain Plasticity. Phy. Rev. Lett. 2006;96:028107. doi: 10.1103/PhysRevLett.96.028107. [PubMed] [Cross Ref]

32. Beggs JM. The criticality hypothesis: how local cortical networks might optimize information processing. Philos. Trans. R. Soc. A. 2008;366:329–343. doi: 10.1098/rsta.2007.2092. [PubMed] [Cross Ref]

33. Moosavi SA, Montakhab A. Mean-field behavior as a result of noisy local dynamics in self-organized criticality: Neuroscience implications. Phys. Rev. E. 2014;89:052139. doi: 10.1103/PhysRevE.89.052139. [PubMed] [Cross Ref]

34. Moosavi SA, Montakhab A. Structural versus dynamical origins of mean-field behavior in a self-organized critical model of neuronal avalanches. Phys. Rev. E. 2015;92:052804. doi: 10.1103/PhysRevE.92.052804. [PubMed] [Cross Ref]

35. de Arcangelis L, Herrmann HJ. Learning as a phenomenon occurring in a critical state. Proc. Natl. Acad. Sci. USA. 2010;107:3977–3981. doi: 10.1073/pnas.0912289107. [PubMed] [Cross Ref]

36. Kinouchi O, Copelli M. Optimal dynamical range of excitable networks at criticality. Nature Physics. 2006;2:384–351. doi: 10.1038/nphys289. [Cross Ref]

37. Larremore DB, Shew WL, Restrepo JG. Predicting criticality and dynamic range in complex networks: effects of topology. Phys. Rev. Lett. 2006;106:058101. doi: 10.1103/PhysRevLett.106.058101. [PubMed] [Cross Ref]

38. Larremore DB, Shew WL, Ott E, Restrepo JG. Effects of network topology, transmission delays, and refractoriness on the response of coupled excitable systems to a stochastic stimulus. Chaos. 2011;21:025117. doi: 10.1063/1.3600760. [PubMed] [Cross Ref]

39. Pei S, et al. How to enhance the dynamic range of excitatory-inhibitory excitable networks. Phys. Rev. E. 2012;86:021909. doi: 10.1103/PhysRevE.86.021909. [PubMed] [Cross Ref]

40. Kandel, E. R., Schwartz, J. H., Jessel, T. M., Siegelbaum, S. A. & Hudspeth, A. J. *Principles of neural science*, 5th ed. (The McGraw-Hill Companies, 2013).

41. Brochini, L. *et al*. Phase transitions and self-organized criticality in networks of stochastic spiking neurons. *Scientific Reports***6**(1) (2016). [PMC free article] [PubMed]

42. Larremore DB, Carpenter MY, Ott E, Restrepo JG. Statistical properties of avalanches in networks. Phys. Rev. E. 2012;85:066131. doi: 10.1103/PhysRevE.85.066131. [PubMed] [Cross Ref]

43. Larremore DB, Shew WL, Ott E, Sorrentino F, Restrepo JG. Inhibition causes ceaseless dynamics in networks of excitable nodes. Phys. Rev. Lett. 2014;112:138103. doi: 10.1103/PhysRevLett.112.138103. [PMC free article] [PubMed] [Cross Ref]

44. Restrepo JG, Ott E, Hunt BR. Approximating the largest eigenvalue of network adjacency matrices. Phys. Rev. E. 2007;76:056119. doi: 10.1103/PhysRevE.76.056119. [PubMed] [Cross Ref]

45. Martin E, Shreim A, Paczuski A. Activity-dependent branching ratios in stocks, solar x-ray flux, and the Bak-Tang-Wiesenfeld sandpile model. Phys. Rev. E. 2010;81:016109. doi: 10.1103/PhysRevE.81.016109. [PubMed] [Cross Ref]

46. Alstrøm P. Mean-field exponents for self-organized critical phenomena. Phys. Rev. A. 1988;38:4905. doi: 10.1103/PhysRevA.38.4905. [PubMed] [Cross Ref]

47. Smith H. On periodic solutions of a delay integral equation modeling epidemics. J. Math. Biol. 1977;4:69–80. doi: 10.1007/BF00276353. [PubMed] [Cross Ref]

48. Hethcote HW, Stech HW, Driessche PVD. Nonlinear oscillations in epidemic models. SIAM J. Appl. Math. 1981;40:1–9. doi: 10.1137/0140001. [Cross Ref]

49. Pinho STR, Andrade RFS. Power-law sensitivity to initial conditions for Abelian directed self-organized critical models. Physica A. 2004;344:601–607. doi: 10.1016/j.physa.2004.06.038. [Cross Ref]

Articles from Scientific Reports are provided here courtesy of **Nature Publishing Group**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |