Collective oscillations of excitable elements: order parameters, bistability and the role of stochasticity
We study the effects of a probabilistic refractory period in the collective behavior of coupled discrete-time excitable cells (SIRS-like cellular automata). Using mean-field analysis and simulations, we show that a synchronized phase with stable collective oscillations exists even with non-deterministic refractory periods. Moreover, further increasing the coupling strength leads to a reentrant transition, where the synchronized phase loses stability. In an intermediate regime, we also observe bistability (and consequently hysteresis) between a synchronized phase and an active but incoherent phase without oscillations. The onset of the oscillations appears in the mean-field equations as a Neimark-Sacker bifurcation, the nature of which (i.e. super- or subcritical) is determined by the first Lyapunov coefficient. This allows us to determine the borders of the oscillating and of the bistable regions. The mean-field prediction thus obtained agrees quantitatively with simulations of complete graphs and, for random graphs, qualitatively predicts the overall structure of the phase diagram. The latter can be obtained from simulations by defining an order parameter suited for detecting collective oscillations of excitable elements. We briefly review other commonly used order parameters and show (via data collapse) that satisfies the expected finite size scaling relations.
- 1 Introduction
- 2 Model
- 3 Complete graph
- 4 Random graph
- 5 Concluding remarks
- A Lyapunov coefficient
Understanding collective oscillations of coupled nonlinear elements remains a challenge from both theoretical and experimental viewpoints. From the theoretical side, much progress has been accomplished since the seminal works of Winfree and Kuramoto [1, 2, 3, 4, 5] on coupled oscillators, upon which recent literature has expanded to include effects of e.g. complex topologies  and noise [7, 8, 9, 10, 11, 12, 13]. From the experimental side, the subject has a longstanding importance in neuroscience: collective neuronal oscillations stood for a long time as candidates for a solution of the so-called binding problem , but emphasis has recently shifted to attentional processes .
Here we are interested in collective oscillations of units which are excitable, i.e. not intrinsically oscillatory. This topic has been experimentally observed in a variety of scenarios, from neuroscience [16, 17] to chemistry , but theoretical approaches have been relatively scarce [19, 20, 21, 22, 23, 24, 25]. In particular, it is not entirely clear to which extent these oscillations are robust with respect to noise. On the one hand, recent studies have shown sufficient conditions for the onset of global oscillations of deterministic excitable units with noisy coupling, emphasizing e.g. the interplay between coupling strength and characteristic time scales of the units [20, 21], or the importance of the topology of the network . On the other hand, the susceptible-infected-recovered-susceptible (SIRS) model on a lattice (i.e. a minimum three-state excitable model) has been thoroughly studied, with several results strongly suggesting that its stochastic Markovian version does not yield sustained oscillations [26, 27, 28, 23] (for non-Markovian models, see e.g. [29, 30]).
This raises the question whether it is possible to find sustained global oscillations in a network of excitable units whose intrinsic dynamics (not only the coupling) is non-deterministic. We therefore propose and study a simple probabilistic model which has a well-defined deterministic limit. To study the phase transitions in the model, we employ an order parameter specifically tailored to assess collective oscillations of excitable systems. These are described in section 2. We study complete graphs as well as random graphs, comparing mean-field calculations with simulations. Results are presented in sections 3 and 4, while section 5 brings our concluding remarks.
2.1 Excitable cellular automata
The minimum model of an excitable system consists of three states, representing quiescence (state 0), excitation (state 1) and refractoriness (state 2)  (a prototypical example being the SIRS model). As shown by Girvan et al., however, such a cyclic three-state deterministic cellular automaton fails to exhibit sustained collective oscillations. For them to become stable, their model needs at least 2 refractory states . To obtain an arbitrary number of refractory states, for each site , let be the consecutive states of the unit (out of the states, the last are refractory , see figure 1).
The cellular automaton version of the probabilistic SIRS model (also called the probabilistic Greenberg-Hastings model ) corresponds to , with intrinsic transitions and governed by constant probabilities, whereas the transition occurs with a probability that usually increases linearly with the number of excited neighbors  (for a study with nonlinear coupling, see ). In the model studied by Girvan et al., on the other hand, all intrinsic transitions (, ) are deterministic.
Here we study an intermediate variant of these models, where all intrinsic transitions are deterministic but the last one, which occurs with probability (see figure 1). The idea is to have a minimum model (lest the number of additional parameters becomes too large) which incorporates non-determinism in the intrinsic dynamics. The choice to make the transition from the last refractory state probabilistic is natural and comes from neuroscience: neuronal dynamics depend on ionic channels which are stochastic [33, 34], so that a neuron may or may not fire when stimulated at the end of its refractory period (the so-called relative refractory period) .
The only transition that still needs to be described is the excitation process . We assume each site is symmetrically connected with other sites. Each active site has a probability of activating a resting neighbour, where is a control parameter (which corresponds to the system branching ratio [35, 22]) and is the average connectivity (). The fraction of active sites
is used to measure network activity at time . At the model shows a transition from an absorbing to an active state, but without sustained oscillations . With deterministic units (), the system undergoes a transition to the oscillatory regime at , which persists indefinitely if is further increased [20, 22].
To motivate the analysis to be developed in section 3, figure 2 shows examples of single-run results in a complete graph with for and increasing values of . For probabilistic units, we observe the transition to an active but nonoscillating state at [Figure 2(a)-(b)] and the second transition to an oscillating state for larger [Figure 2(b)-(c)]. Contrary to what is observed in the deterministic model, however, in the probabilistic model this second transition is reentrant with respect to the coupling strength , as shown in figure 2(c)-(d).
The nature of this reentrant transition will be clarified in section 3. In order to analyze it properly, though, one needs to define an adequate order parameter to detect synchronization among excitable elements.
2.3 Order parameters
where . Note that corresponds to the center of mass of the phases of the units. This works fine when the system is composed of coupled uniform oscillators, because rotational symmetry will ensure that, in the absence of sustained collective oscillations, the order parameter vanishes in the thermodynamic limit. Consider, however, the present case of excitable elements. The trivial absorbing state (, ), which is always a collective solution of the dynamics, yields a nonzero (in fact, maximum!) Kuramoto order parameter. Indeed, in the absorbing state units are all “perfectly synchronized” in the sense that they always have the same state. But this is clearly not what one wants to detect. This problem persists for (below the onset of collective oscillations), where a small fraction of the units are active (on average), whereas most remain quiescent. In that case, has a constant bias towards the absorbing state, around which it will fluctuate.
Collective oscillations correspond to rotations of which may be misdetected in the averaging procedure owing to a lurking constant vector. One possibility which has been used to avoid the weight of the absorbing state is excluding terms with from the sum in eq. (2) [19, 21]. Another strategy makes use of the standard deviation (measured along time) of to detect oscillations , though neither procedure is easily extensible to continuous-phase systems.
To account for a system of continuous-phase units which may have an arbitrary number of preferred phases, one could employ the angular momentum , where [37, 38]. In our cellular automata, this would require the discretization of the time derivative, which could introduce unnecessary numerical errors. Alternatively, Shinomoto and Kuramoto have previously proposed
which amounts to subtracting the constant bias from before the averaging procedure . However, this can be computationally expensive, requiring the storage of the whole time series. Making use of a similar idea, but with much less computational bookkeeping, in the following we will characterize collective oscillations via the order parameter
which can be seen as a generalized standard deviation of . Differently from , obtaining is computationally inexpensive, as the means over time are now separated and may be calculated along with the simulation. In section 4.3 we will show that satisfies scaling relations near a phase transition, as expected for an order parameter.
One may grasp intuition about by considering the different time series in figure 2. Let be the stationary value , which is an order parameter in its own right, measuring whether or not the network is active . In figs. 2(a) and (b), although the system goes from to as the system goes from an absorbing to an active fixed point, both have , because there are no oscillations after the transient. In figure 2(c), on the other hand, we have both and , as the oscillatory state becomes stable after increasing . Finally, for figure 2(d), the oscillations are unstable, and .
3 Complete graph
We start with the complete graph because it is presumably the topology most prone to exhibiting stable collective oscillations. Besides, it allows a comparison between simulations and analytical results (see below). From now on, we will focus on the effect of the probabilistic dynamics () on synchronization and will fix .
We have simulated complete graphs () with sizes varying from to . At each time step, intrinsic transitions occur as described in section 2. The transition is governed by the number of active () sites at time , each of which can activate a quiescent cell with probability . Therefore, the probability of a quiescent cell being activated by at least one of its active neighbours is
which renders the simulations relatively simple despite the large system sizes. Finally, initial conditions must be chosen to avoid large amplitudes during the transient, which could throw the system into the absorbing state . Apart from that, the effects we report below are robust with respect to the initial conditions and we have arbitrarily fixed and .
To illustrate the kind of reentrant transition exemplified in figure 2, we show in figure 3(a) the order parameter as a function of the coupling parameter . Starting at some , for each value of we let the system evolve during a transient of time steps, after which we start measuring the order parameter up to time steps. We then increase by a constant amount and repeat the procedure, with the initial condition of the system for each value corresponding to the final condition of the preceding value. Both constants are chosen so the rate of change is small (). After a maximum value is reached, is sequentially decreased by the same amount down to .
As shown in figure 3(a), the reentrance of the transition to collective oscillations is captured by the non-monotonic behavior of , which departs from zero at some lower value and returns to zero at some upper value . Moreover, we have found that while the first transition is always continuous, the second transition can be discontinuous. The fingerprint of the discontinuity is the hysteresis observed in the order parameter: above , the only stable state of the system has constant (but nonzero) , thus no oscillations (). If we decrease , oscillations do not reappear at , but rather at a lower value . There is therefore a region of bistability in parameter space where collective oscillations () can coexist with an active () but non-oscillating () state. As it turns out [see Figure 3(a)], for the size of this bistable region is rather sensitive to the system size. Smaller systems tend to be perturbed away from collective oscillations by larger fluctuations, leading to smaller hysteresis cycles.
As is decreased, the mean and variance of the refractory periods of the units increase, rendering the whole system noisier. This might be the explanation for the result in figure 3(b), which shows a smaller reentrant region with oscillations (). The width of the hysteresis cycle also decreases, with both and decreasing. Furthermore, (along with the width of the hysteresis cycle) becomes less sensitive to the system size, which could be due to the variance of the refractory periods overcoming the effects of small-size fluctuations. Albeit subtly, also increases slowly with decreasing , as will be seen in figure 4.
Further decreasing [Figure 3(c)], the bistable region vanishes, whereas the transition to a collectively oscillating state remains. Finally, for sufficiently small , collective oscillations are no longer stable [Figure 3(d)].
In the following, we will show that these transitions can be quantitatively reproduced by a low-dimensional mean-field analysis. For a controlled comparison, and were heuristically defined as the values of (averaged over runs) where first rose above some threshold value , respectively for increasing and decreasing values of . On the other hand, was defined as the value of (increasing) where first fell below .
3.2 Mean-field analysis
In the following we apply the standard mean-field (MF) approximation  to the equations governing our system. We follow closely the steps of Refs. [20, 41, 22], where every site is considered to have neighbors, a fraction of which is excited at time . If a given site is at rest [with probability ], the probability of it becoming excited by at least one of its excited neighbors is
The dynamics of the system is then described by the following closed set of equations:
where the normalization condition
renders (7) redundant and reduces the system to a -dimensional map . Therefore, increasing the duration of the refractory period amounts to an increase in the complexity of the mean-field calculations.
For the complete graph, and mean field is exact. In the thermodynamic limit, (6) becomes
In other words, in the mean field approach the state of the system is completely described by , which evolves according to . Note that is the only component of which is nonlinear [see (13)].
which can be numerically solved. Expanding near (note that is always a solution), one easily obtains the transition from an absorbing to an active (steady) state at . The collective oscillations appear when the active state becomes unstable.
3.3 Linear stability
Considering a small perturbation such that , the linearized dynamics can be written as , where is the Jacobian matrix calculated at the fixed point:
and . The eigenvalues of determine whether is stable () or unstable () (for simplicity, in the following we employ , where ).
We expect to pinpoint the transition to the oscillatory state by looking for a Neimark-Sacker (NS) bifurcation in the mean field equations (which is the discrete-time analog of the Andronov-Hopf (AH) bifurcation in continuous time ). In other words, for fixed , we have with at . The relation between and the pair of critical values depicted in figure 3(a) will be clarified below.
Like the AH bifurcation, the NS bifurcation also comes in two different flavours: in the supercritical case, a stable closed invariant curve (CIC — the discrete-time analog of a limit cycle) is born at and grows continually from zero amplitude; in the subcritical case, an unstable CIC exists below and engulfs the fixed point at (above which the system is typically attracted to another pre-existing, but stable, CIC). Since the order parameter increases with the amplitude of the oscillations (which is, roughly speaking, proportional to ), a supercritical (subcritical) NS bifurcation in the mean-field equations is suggestive of a continuous (discontinuous) phase transition in the system (see e.g. ).
The sign of the first Lyapunov coefficient  indicates if the Neimark-Sacker bifurcation is supercritical () or subcritical (). Its calculation is briefly reviewed in the Appendix. We now have the necessary tools to unveil the complete phase diagram.
3.4 Phase diagram
We have run simulations of complete graphs with excitable units and employed the protocol described in section 3.1 with to detect the width of the hysteresis loop (coexistence region). We have tested and verified that longer values of do not change our results significantly. The phase diagram thus obtained from the simulations is shown with symbols in figure 4 (the horizontal gray lines show the values of used in figure 3). To obtain the NS lines of the mean-field equations, we have numerically explored a special test function  , which changes its sign at the NS bifurcation. As the first Lyapunov coefficient determines whether the bifurcation is super- or subcritical, its value along the bifurcation line is shown in the inset of figure 4 (changing sign at ). The solid lines in figure 4 show the supercritical (red) and subcritical (blue) bifurcation curves where .
Comparing the solid lines with the symbols in figure 4, we observe that linear stability analysis accurately accounts for the transition from an active (but non-oscillating) phase to an oscillating phase. In other words, it correctly predicts the lines and , where in both cases the non-oscillating phase loses stability.
The transition at , however, cannot be predicted by linear analysis. Note that in this case it is the stable CIC that loses its stability (that of the fixed point remaining intact). This hints at the existence of a global bifurcation, which can be numerically detected in the MF equations by direct iteration of the map determined by equations (9), (10) and (13). To compare the -dimensional MF map with system simulations, we rewrite the complex vector from (2) and (11) as
where . We can thus calculate for the MF map and subject it to the same protocol used for detecting the coexistence region in the simulations. The black solid line in figure 4 shows the obtained by iteration of the map, which is in good agreement with simulations (symbols).
Note that in the lower part of the oscillating phase () the order parameter detects oscillations in the simulations which are not predicted by the mean-field analysis. This phenomenon is due to stochastic oscillations, as recently explained by Risau-Gusman and Abramson : the fixed point in the conflicting region is in fact stable, but with an eigenvalue with a nonzero imaginary part. Inevitable fluctuations throw the system away from the stable point, to which it returns in spiral-like trajectories, yielding a nonzero even for very large system sizes [43, 23].
For , we recover a quenched variant of the model by Girvan et al. . In this regime where intrinsic transitions are deterministic, increasing the coupling will only reinforce collective oscillations, and the fixed point never regains stability (i.e. ). This suggests that even small amounts of noise in the intrinsic dynamics of excitable elements can lead to qualitatively different collective behavior in a regime of strong coupling.
4 Random graph
4.1 Mean-field and simulation results
To understand network topology effects on synchronization, we study a bidirectional random graph similar to an Erdős-Rényi  network, where links connect randomly chosen pairs  and remain frozen (“quenched”) throughout each run (and in each run, a new realization of the network is created). The main difference with respect to the complete graph (CG) is that in the random graph (RG) the value of is bounded from above: [see (6)]. The mean-field calculations, however, are otherwise similar to that of the complete graph, with (6) replacing (12). We therefore applied to the RG problem the same procedures for determining the stability of the solutions, the nature of the NS bifurcation and the boundary of the bistability region (see section 3).
Given their uncorrelated assigment of links and short distances among sites, random graphs are usually regarded as the natural topology in which mean-field predictions are expected to hold. Indeed, simulations and mean-field calculations agree nearly perfectly as far as the phase transition at  is concerned. In figure 5(a)-(b) we see that a good agreement is also observed in the transitions to a synchronized phase for large values of . Note, however, that simulations and mean-field predictions differ at the rightmost boundary of the bistable region, and the disagreement worsens as decreases. As shown in figure 5(c), for smaller values of bistability was not even detected in the simulations, and the oscillating phase is substantially smaller than predicted by mean field.
4.2 Annealed random graphs
Could correlations (which are neglected by the mean-field approximation) account for the discrepancy observed in figure 5? In order to assess the role of the correlations associated with the quenched connectivity, we studied an annealed variant of the model where the neighbors of each site are randomly chosen at each time step [20, 45]. Results are shown in figure 6, in which we restrict ourselves to smaller values of because results are essentially indistinguishable from the quenched case for large .
For smaller , three features are noteworthy. First, the agreement with mean field results is recovered (apart from the stochastic oscillations in the lower end of the oscillating phase, like in the previous cases). This therefore confirms the suspicion that correlations associated to the quenched connectivity can indeed undermine collective oscillations. This is not surprising, given the difficulty of establishing collective oscillations of excitable elements in hypercubic lattices [19, 23].
Second, the bound impoverishes the repertoire of detected phenomena for small [see the forbidden gray regions in figure 6(b) and (c)]. Note that the coexistence region shrinks as decreases. In fact, since varies very slowly with [see the inset of figure 6(c)], there is a minimum value of [satisfying ] below which the system shows no bistability (i.e. there is no change in the sign of ). We have numerically estimated .
Finally, note that the oscillating phase extends into lower values of for decreasing . This is a rather counter-intuitive result. It means that, for fixed and , it is possible to take the system from a non-oscillating to an oscillating phase by lowering the connectivity . This is a particularity of the annealed RG (and the mean-field approximation), however. Note in figure 5 that the opposite (and expected) trend is observed for quenched random graphs. It remains to be studied whether refining the approximation (including e.g. first-neighbor correlations [41, 28]) can reconcile mean-field with quenched random graph results.
4.3 Finite-size scaling
Our phase diagrams rely heavily on the proposed order parameter . The inset of figure 7 indicates that near its behavior becomes increasingly abrupt with increasing system size . In order to confirm that indeed possesses the basic properties of a bona fide order parameter, here we show that at the supercritical NS bifurcation (i.e. a second order phase transition) it satisfies the scaling relations that would be expected from standard finite-size scaling (FSS) theory.
Defining , FSS predicts that for a lattice with linear size , where the critical exponents are defined as and in the limit (where is the correlation length). This holds for below the upper critical dimension . For , mean-field exponents are expected and the scaling relation has to be modified  with , so
This modified version of the usual FSS relation was shown to also hold for infinitely coordinated networks  (i.e. the complete graph).
For large argument, the scaling function in Equation (17) becomes , as usual . In the subcritical regime (), one expects , so . For this to be true, when . Figure 7 shows an excellent data collapse for (quenched) random graphs with different system sizes. Consistent power laws are obtained with standard mean-field exponents (as expected for random graphs), namely , and . Similar results are obtained with complete graphs (not shown).
5 Concluding remarks
We have studied the effects of a probabilistic refractory period in the collective behavior of a large number of coupled excitable cellular automata. We have obtained the mean-field solution of the model and compared it with simulations of complete as well as random graphs. The continuous phase transition to a synchronized regime is associated to a Neimark-Sacker bifurcation in the mean-field equations.
This scenario is similar to what has been previously obtained by Girvan et al.  in a model of deterministic excitable automata ( in our model). The effects of setting , however, are drastic, and appear for sufficiently strong coupling , when we have observed that oscillations vanish. This is in contrast with the transition to the absorbing state found in Ref.  for strong coupling. While in their model the transition is due to very large amplitudes driving the system into rest (), in our model the system is thrown into an active albeit disordered state, which still has , but no oscillations.
Furthermore, only for non-deterministic excitable elements do we observe bistability, with an oscillating and an active (but non-oscillating) phase coexisting. This leads to hysteresis cycles, whose sizes can depend on the system size (notably for the complete graph) and eventually disappear in random graphs with small enough mean connectivity .
Although we have restricted ourselves to , preliminary results suggest that the overall scenario is preserved for larger values of , specially regarding the first transition at . The observation of bistability and hysteresis is more difficult for larger values of , owing to stronger finite-size effects. These are similar to those reported by Girvan et al.: for finite and sufficiently strong coupling, the CIC grows in amplitude and nears the absorbing state, to which the system is thrown by fluctuations . Distinguishing between that type of transition and the bistability reported here is not obvious and remains to be studied.
It is interesting to note that a simple model allows the straightforward application of standard techniques of nonlinear dynamics to the study of these phase transitions, which could otherwise be difficult to tackle. Note that in the limit , each excitable unit in our model has, at the end of its refractory period, a perfect memory of its past time steps. Incorporating this memory in a continuous-time model would require non-Markovian dynamics, as recently proposed by Gonçalves et al. . Interestingly, their model also shows an oscillating phase with reentrance, whose size decreases as memory time decreases. In our model, this corresponds to lowering , with a similar effect on the collective behavior. It remains to be investigated whether bistability also occurs in non-Markovian continuous-time models.
Finally, a note of caution is in order regarding the scaling results of figure 7. The fact that a data collapse is obtained employing an upper critical dimension by no means implies that a lower critical dimension exists. In fact, we are not aware of models which exhibit collective oscillations of excitable elements in hypercubic lattices (though they do appear if stimulated with a Poisson drive ). It remains to be investigated whether the results of this model hold for small-world networks and other complex topologies [19, 21], which are extremely appealing for applications in Neuroscience. This is currently under investigations (results will be published elsewhere).
In summary, we have shown that collective oscillations of excitable elements have robustness to a certain degree of stochasticity in their intrinsic dynamics. For fixed coupling, there is a critical value of , below which no oscillations are stable. Fixing and increasing the coupling , on the other hand, leads to interesting new phenomena such as bistability and discontinuous transitions. Taken together, our results suggest that even weakly noisy dynamics can qualitatively change the collective oscillating behavior. Studies attempting to verify whether these phenomena are observed in more detailed excitable networks (e.g. modelled by stochastic differential equations) would certainly be welcome.
Appendix A Lyapunov coefficient
Let and be respectively the right and left (adjoint) eigenvector of the Jacobian matrix:
with both normalized: (brackets denote the standard complex inner product). Let also and be multilinear functions proportional to the first nonlinear terms of the Taylor expansion of at , i.e.,
If we now define
where is the identity matrix and is the conjugate of , the coefficient is finally given by 
-  A. T. Winfree. Biological rhythms and the behavior of populations of coupled oscillators. J. Theor. Biol., 16(1):15–42, 1967.
-  Y. Kuramoto, editor. Chemical Oscillations, Waves and Turbulence. Dover, Berlin, 1984.
-  S. H. Strogatz. From Kuramoto to Crawford: Exploring the onset of synchronization in populations of coupled oscillators. Physica D, 143:1–20, 2000.
-  A. Pikovsky, M. Rosenblum, and J. Kurths. Synchronization: A Universal Concept in Nonlinear Sciences. Cambridge University Press, Cambridge, UK, 2001.
-  J. Acebrón, L. Bonilla, C. Pérez Vicente, F. Ritort, and R. Spigler. The Kuramoto model: A simple paradigm for synchronization phenomena. Rev. Mod. Phys., 77(1):137–185, 2005.
-  A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou. Synchronization in complex networks. Phys. Rep., 469(3):93–153, 2008.
-  T. Risler, J. Prost, and F. Jülicher. Universal critical behavior of noisy coupled oscillators. Phys. Rev. Lett., 93(17):175702, 2004.
-  T. Risler, J. Prost, and F. Jülicher. Universal critical behavior of noisy coupled oscillators: A renormalization group study. Phys. Rev. E, 72(1):016130, 2005.
-  K. Wood, C. Van den Broeck, R. Kawai, and K. Lindenberg. Universality of synchrony: Critical behavior in a discrete model of stochastic phase-coupled oscillators. Phys. Rev. Lett., 96:145701, 2006.
-  K. Wood, C. Van den Broeck, R. Kawai, and K. Lindenberg. Critical behavior and synchronization of discrete stochastic phase-coupled oscillators. Phys. Rev. E, 74(3):031113, 2006.
-  K. Wood, C. Van den Broeck, R. Kawai, and K. Lindenberg. Effects of disorder on synchronization of discrete phase-coupled oscillators. Phys. Rev. E, 75(4):041107, 2007.
-  K. Wood, C. Van den Broeck, R. Kawai, and K. Lindenberg. Continuous and discontinuous phase transitions and partial synchronization in stochastic three-state oscillators. Phys. Rev. E, 76(4):041132, 2007.
-  H. Khoshbakht, F. Shahbazi, and K. A. Samani. Phase synchronization on scale-free and random networks in the presence of noise. J. Stat. Mech., 2008(10):P10020, 2008.
-  W. Singer. Neuronal synchrony: a versatile code for the definition of relations? Neuron, 24(1):49–65, 1999.
-  P. J. Uhlhaas, G. Pipa, B. Lima, L. Melloni, S. Neuenschwander, D. Nikolić, and W. Singer. Neural synchrony in cortical networks: history, concept and current status. Front. Integr. Neurosci., 3:17, 2009.
-  W. Shew, H. Yang, T. Petermann, R. Roy, and D. Plenz. Neuronal avalanches imply maximum dynamic range in cortical networks at criticality. J. Neurosci., 29(49):15595–15600, 2009.
-  C. H. Ko, Y. R. Yamada, D. K. Welsh, E. D. Buhr, A. C. Liu, E. E. Zhang, M. R. Ralph, S. A. Kay, D. B. Forger, and J. S. Takahashi. Emergence of Noise-Induced Oscillations in the Central Circadian Pacemaker. PLoS Biology, 8(10):e1000513, 2010.
-  A. F. Taylor, M. R. Tinsley, F. Wang, Z. Huang, and K. Showalter. Dynamical quorum sensing and synchronization in large populations of chemical oscillators. Science, 323(5914):614–7, 2009.
-  M. Kuperman and G. Abramson. Small world effect in an epidemiological model. Phys. Rev. Lett., 86(13):2909–2912, 2001.
-  M. Girvan, D. S. Callaway, M. E. J. Newman, and S. H. Strogatz. Simple model of epidemics with pathogen mutation. Phys. Rev. E, 65:031915, 2002.
-  P. Gade and S. Sinha. Dynamic transitions in small world networks: Approach to equilibrium limit. Phys. Rev. E, 72(5):1–4, 2005.
-  O. Kinouchi and M. Copelli. Optimal dynamical range of excitable networks at criticality. Nat. Phys., 2:348–351, 2006.
-  V. R. V. Assis and M. Copelli. Discontinuous nonequilibrium phase transitions in a nonlinearly pulse-coupled excitable lattice model. Phys. Rev. E, 80:061105, 2009.
-  Q.-X. Liu, R.-H. Wang, and Z. Jin. Persistence, extinction and spatio-temporal synchronization of SIRS spatial models. J. Stat. Mech., 2009(07):P07007, 2009.
-  P. McGraw and M. Menzinger. Self-Sustaining Oscillations in Complex Networks of Excitable Elements. arXiv:1009.0537v1 [cond-mat.dis-nn], 2010.
-  G. Rozhnova and A. Nunes. Fluctuations and oscillations in a simple epidemic model. Phys. Rev. E, 79:041922, 2009.
-  G. Rozhnova and A. Nunes. SIRS dynamics on random networks: Simulations and analytical models. In J. Zhou, editor, Complex Sciences, volume 4, pages 792–797. Springer Berlin Heidelberg, 2009. arXiv:0812.1812v1 [q-bio.PE].
-  G. Rozhnova and A. Nunes. Cluster approximations for infection dynamics on random networks. Phys. Rev. E, 80:051915, 2009.
-  B. Lindner, J. García-Ojalvo, A. Neiman, and L. Schimansky-Geier. Effects of noise in excitable systems. Phys. Rep., 392:321–424, 2004.
-  S. Gonçalves, G Abramson, and M. F. C. Gomes. Oscillations in SIRS model with distributed delays. arXiv:0912.1250v3 [q-bio.PE], 2009.
-  J. M. Greenberg and S. P. Hastings. Spatial patterns for discrete models of diffusion in excitable media. SIAM J. Appl. Math., 34:515–523, 1978.
-  V. R. V. Assis and M. Copelli. Dynamic range of hypercubic stochastic excitable media. Phys. Rev. E, 77:011923, 2008.
-  C. Koch. Biophysics of Computation. Oxford University Press, New York, 1999.
-  P. V. Carelli, M. B. Reyes, J. C. Sartorelli, and R. D. Pinto. Whole cell stochastic model reproduces the irregularities found in the membrane potential of bursting neurons. J. Neurophysiol., 94:1169–1179, 2005.
-  T. E. Harris. The Theory of Branching Processes. Springer, 1963.
-  J. Garcia-Ojalvo, M. B. Elowitz, and S. H. Strogatz. Modeling a synthetic multicellular clock: repressilators coupled by quorum sensing. Proc. Natl. Acad. Sci. USA, 101(30):10955–60, 2004.
-  Y. Kuramoto, T. Aoyagi, I. Nishikawa, T. Chawanya, and K. Okuda. Neural Network Model Carrying Phase Information with Application to Collective Dynamics. Prog. Theor. Phys., 87:1119–1126, 1992.
-  H. Ohta and S. Sasa. Critical phenomena in globally coupled excitable elements. Phys. Rev. E, 78(6):065101, 2008.
-  S. Shinomoto and Y. Kuramoto. Phase transitions in active rotator systems. Progr. Theoret. Phys., 75(5):1105–1110, 1986.
-  J. Marro and R. Dickman. Nonequilibrium Phase Transition in Lattice Models. Cambridge University Press, Cambridge, 1999.
-  L. S. Furtado and M. Copelli. Response of electrically coupled spiking neurons: a cellular automaton approach. Phys. Rev. E, 73:011907, 2006.
-  Y.A. Kuznetsov. Elements of Applied Bifurcation Theory (Applied Mathematical Sciences). Springer, 2nd edition, 1998.
-  S. Risau-Gusman and G. Abramson. Bounding the quality of stochastic oscillations in populations models. Eur. Phys. J. B, 60:515–520, 2007.
-  P. Erdős and A. Rényi. On the evolution of random graphs. Publ. Math. Inst. Hung. Acad. Sci, 5:17–61, 1960.
-  S. Sinha, J. Saramäki, and K. Kaski. Emergence of self-sustained patterns in small-world excitable media. Phys. Rev. E, 76(1):015101, July 2007.
-  E. Brézin. An investigation of finite size scaling. J. Phys.(Paris), 43(1):15–22, 1982.
-  R. Botet, R. Jullien, and P. Pfeuty. Size scaling for infinitely coordinated systems. Phys. Rev. Lett., 49(7):478–481, 1982.
-  T. J. Lewis and J. Rinzel. Self-organized synchronous oscillations in a network of excitable cells coupled by gap junctions. Network: Comput. Neural Syst., 11:299–320, 2000.