{'fraction_non_alphanumeric': 0.044757072103552, 'fraction_numerical': 0.014280939503509464, 'mean_word_length': 4.42296918767507, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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 q suited for detecting collective oscillations of excitable elements. We briefly review other commonly used order parameters and show (via data collapse) that q satisfies the expected finite size scaling relations.


Introduction
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 [6] 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 [14], but emphasis has recently shifted to attentional processes [15].
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 [18], 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 [19]. On the other hand, the susceptibleinfected-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.

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) [29] (a prototypical example being the SIRS model). As shown by Girvan et al., however, such a cyclic threestate deterministic cellular automaton fails to exhibit sustained collective oscillations. For them to become stable, their model needs at least 2 refractory states [20]. To obtain an arbitrary number of refractory states, for each site j = 1, ..., N, let s j = 0, 1, 2, ..., τ be the consecutive states of the unit (out of the τ +1 states, the last τ −1 are refractory [20], see figure 1).
The cellular automaton version of the probabilistic SIRS model (also called the probabilistic Greenberg-Hastings model [31]) corresponds to τ = 2, with intrinsic transitions 1 → 2 and 2 → 0 governed by constant probabilities, whereas the transition 0 → 1 occurs with a probability that usually increases linearly with the number of excited neighbors [32] (for a study with nonlinear coupling, see [23]). In the model studied by Girvan et al., on the other hand, all intrinsic transitions (s j → (s j + 1) mod (τ + 1), s j = 0) 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 p γ (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) [33]. Single-cell dynamics. p inf is the probability of activation (0 → 1) from neighbours, p γ is the probability of transitioning from the relative refractory state (τ ) to the rest state (0). Light gray states are refractory. All other transitions are deterministic.

Coupling
The only transition that still needs to be described is the excitation process 0 → 1. We assume each site j is symmetrically connected with k j other sites. Each active site has a probability σ/K of activating a resting neighbour, where σ is a control parameter (which corresponds to the system branching ratio [35,22]) and K is the average connectivity (K = k j ). The fraction of active sites is used to measure network activity at time t. At σ = 1 the model shows a transition from an absorbing to an active state, but without sustained oscillations [22]. With deterministic units (p γ = 1), the system undergoes a transition to the oscillatory regime at σ = σ c (p γ = 1) > 1, 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 N = 5 × 10 5 for p γ = 0.85 and increasing values of σ. For probabilistic units, we observe the transition to an active but nonoscillating state at σ = 1 [ Figure 2 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.

Order parameters
Most studies of synchronization employ the Kuramoto order parameter [2,3], which corresponds to the time and ensemble average of the norm of the complex vector where θ j = 2πs j /(τ + 1). Note that Z 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 (s j = 0, ∀j), 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 σ 1 (below the onset of collective oscillations), where a small fraction of the units are active (on average), whereas most remain quiescent. In that case, Z has a constant bias towards the absorbing state, around which it will fluctuate.
Collective oscillations correspond to rotations of Z 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 s j = 0 from the sum in eq. (2) [19,21]. Another strategy makes use of the standard deviation (measured along time) of P t (1) to detect oscillations [36], 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 L ≡ X∂ t Y − Y ∂ t X, where Z = X + iY [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 Z before the averaging procedure [39]. 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 Z(t). Differently fromq, obtaining q 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 q satisfies scaling relations near a phase transition, as expected for an order parameter.
One may grasp intuition about q by considering the different time series in figure 2. Let P * 1 be the stationary value lim t→∞ P t (1) t , which is an order parameter in its own right, measuring whether or not the network is active [22]. In figs. 2(a) and (b), although the system goes from P * 1 = 0 to P * 1 = 0 as the system goes from an absorbing to an active fixed point, both have q = 0, because there are no oscillations after the transient. In figure 2(c), on the other hand, we have both P * 1 = 0 and q > 0, as the oscillatory state becomes stable after increasing σ. Finally, for figure 2(d), the oscillations are unstable, P * 1 = 0 and q = 0.

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 (p γ ) on synchronization and will fix τ = 3.

Simulations
We have simulated complete graphs (k j = K = N − 1) with sizes varying from N = 10 5 to N = 10 6 . At each time step, intrinsic transitions occur as described in section 2. The transition 0 → 1 is governed by the number N 1 (t) of active (s i = 1) sites at time t, each of which can activate a quiescent cell with probability σ/N. Therefore, the probability of a quiescent cell being activated by at least one of its N 1 (t) 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 [20]. Apart from that, the effects we report below are robust with respect to the initial conditions and we have arbitrarily fixed P 0 (0) = 0.8 and P 0 (1) = 0.2.
To illustrate the kind of reentrant transition exemplified in figure 2, we show in figure 3(a) the order parameter q as a function of the coupling parameter σ. Starting at some σ min > 1, for each value of σ we let the system evolve during a transient of t trans time steps, after which we start measuring the order parameter q up to t = t max 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 (δσ/t max ≪ 1). After a maximum value σ max is reached, σ is sequentially decreased by the same amount δσ down to σ min .
As shown in figure 3(a), the reentrance of the transition to collective oscillations is captured by the non-monotonic behavior of q(σ), which departs from zero at some lower value σ c (p γ ) and returns to zero at some upper value σ c 2 (p γ ). 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 σ c 2 , the only stable state of the system has constant (but nonzero) P t (1), thus no oscillations (q = 0). If we decrease σ, oscillations do not reappear at σ c 2 , but rather at a lower value σ c 1 . There is therefore a region of bistability σ ∈ [σ c 1 , σ c 2 ] in parameter space where collective oscillations (q > 0) can coexist with an active (P t (1) > 0) but non-oscillating (q = 0) state. As it turns out [see Figure 3(a)], for p γ = 0.95 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 p γ 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 (q > 0). The width of the hysteresis cycle also decreases, with both σ c 1 and σ c 2 decreasing. Furthermore, σ c 2 (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, σ c also increases slowly with decreasing p γ , as will be seen in figure 4.
Further decreasing p γ [ Figure 3(c)], the bistable region vanishes, whereas the transition to a collectively oscillating state remains. Finally, for sufficiently small p γ , 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, σ c and σ c 1 were heuristically defined as the values of σ (averaged over n runs) where q first rose above some threshold value q min ∝ 1/ √ N, respectively for increasing and decreasing values of σ. On the other hand, σ c 2 was defined as the value of (increasing) σ where q first fell below q min .

Mean-field analysis
In the following we apply the standard mean-field (MF) approximation [40] to the equations governing our system. We follow closely the steps of Refs. [20,41,22], where every site is considered to have K neighbors, a fraction P t (1) of which is excited at time t. If a given site is at rest [with probability P t (0)], 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 [41]. Therefore, increasing the duration of the refractory period amounts to an increase in the complexity of the mean-field calculations.
For the complete graph, K = N − 1 and mean field is exact. In the thermodynamic limit, (6) becomes lim N →∞ which, from (7) and (11), leads to In other words, in the mean field approach the state of the system is completely described by P t ≡ (P t (1), P t (2), . . . , P t (τ )) T , which evolves according to P t+1 = F ( P t ). Note that F 1 is the only component of F which is nonlinear [see (13)]. From (9), the fixed point for any 2 ≤ s ≤ τ −1 is clearly P * s ≡ P ∞ (s) = P ∞ (s−1) = · · · = P * 1 which, upon substitution in (10), gives P * τ = P * 1 /p γ . Finally, in its steady state, (13) becomes: which can be numerically solved. Expanding near P * 1 ≃ 0 (note that P * 1 = 0 is always a solution), one easily obtains the transition from an absorbing to an active (steady) state at σ = 1 [22]. The collective oscillations appear when the active state becomes unstable.

Linear stability
Considering a small perturbation η t (s) such that P t (s) = P * s + η t (s), the linearized dynamics can be written as η t+1 = A η t , where A ij = ∂F i /∂P t (j)| P * is the Jacobian matrix calculated at the fixed point: and g(σ, P * 1 ) = e −σP * of A determine whether P * is stable (max j |µ j | < 1) or unstable (max j |µ j | > 1) (for simplicity, in the following we employ µ ≡ µ k , where k = argmax j |µ j |). 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 discretetime analog of the Andronov-Hopf (AH) bifurcation in continuous time [42]). In other words, for fixed p γ , we have |µ| = 1 with Im(µ) = 0 at σ = σ N S . The relation between σ N S and the pair {σ c 1 , σ c 2 } 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 σ N S and grows continually from zero amplitude; in the subcritical case, an unstable CIC exists below σ N S and engulfs the fixed point P * at σ N S (above which the system is typically attracted to another pre-existing, but stable, CIC). Since the order parameter q 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. [23]).
The sign of the first Lyapunov coefficient l 1 [42] indicates if the Neimark-Sacker bifurcation is supercritical (l 1 < 0) or subcritical (l 1 > 0). Its calculation is briefly reviewed in the Appendix. We now have the necessary tools to unveil the complete phase diagram.

Phase diagram
We have run simulations of complete graphs with N = 10 6 excitable units and employed the protocol described in section 3.1 with t max = 10 3 to detect the width of the hysteresis loop (coexistence region). We have tested and verified that longer values of t max 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 p γ used in figure 3). To obtain the NS lines of the mean-field equations, we have numerically explored a special test function [42] Φ N S (σ; p γ ) = τ m<n (1 − µ n µ m ), which changes its sign at the NS bifurcation. As the first Lyapunov coefficient l 1 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 σ T ). The solid lines in figure 4 show the supercritical (red) and subcritical (blue) bifurcation curves where Φ N S (σ N S ; p γ ) = 0.
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 nonoscillating) phase to an oscillating phase. In other words, it correctly predicts the lines σ c (p γ ) and σ c 1 (p γ ), where in both cases the non-oscillating phase loses stability. The transition at σ c 2 (p γ ), 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 P * remaining intact). This hints at the existence of a global bifurcation, which can be  . Phase diagram for the complete graph (triangles) with N = 10 6 (mean over 50 runs, t max = 10 3 , t trans = 5 × 10 2 ). Solid red (blue) line is the supercritical (subcritical) Neimark-Sacker bifurcation predicted from linear analysis. Black line marks the (discontinuous) stability limit of the oscillating phase (as predicted by mean field). Note that the black and blue lines approach each other very closely before merging at σ T . Inset: first Lyapunov coefficient l 1 (σ).
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 Z from (2) and (11) as where φ s = 2πs/(τ + 1). We can thus calculate q 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 σ c 2 (p γ ) obtained by iteration of the map, which is in good agreement with simulations (symbols).
Note that in the lower part of the oscillating phase (p γ 0.8) 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 [43]: 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 q even for very large system sizes [43,23].
For p γ = 1, we recover a quenched variant of the model by Girvan et al. [20]. In this regime where intrinsic transitions are deterministic, increasing the coupling will only reinforce collective oscillations, and the fixed point P * 1 never regains stability (i.e. σ c 1 → ∞). 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.

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 [44] network, where NK/2 links connect randomly chosen pairs [22] 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: σ ≤ K [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 σ = 1 [22] 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 K. Note, however, that simulations and mean-field predictions differ at the rightmost boundary of the bistable region, and the disagreement worsens as K decreases. As shown in figure 5(c), for smaller values of K bistability was not even detected in the simulations, and the oscillating phase is substantially smaller than predicted by mean field.

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 K 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 K because results are essentially indistinguishable from the quenched case for large K.
For smaller K, 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 Neimark-Sacker bifurcation predicted from linear analysis. Black line marks the (discontinuous) stability limit of the oscillating phase (as predicted by mean field). Symbols are obtained from simulations (mean over 5 runs) with t max = 3 × 10 3 steps (t trans = 2 × 10 3 steps) and N = 10 5 . Standard errors are smaller than symbol size. The purple dashed line is a guide for the eyes and marks the stability limit of the oscillating phase for simulations.
oscillations. This is not surprising, given the difficulty of establishing collective oscillations of excitable elements in hypercubic lattices [19,23].
Second, the bound σ < K impoverishes the repertoire of detected phenomena for small K [see the forbidden gray regions in figure 6(b) and (c)]. Note that the coexistence region shrinks as K decreases. In fact, since σ T varies very slowly with K [see the inset of figure 6(c)], there is a minimum value of K [satisfying K c = σ T (K c )] below which the system shows no bistability (i.e. there is no change in the sign of l 1 ). We have numerically estimated K c = 7.8074 (1).
Finally, note that the oscillating phase extends into lower values of p γ for decreasing K. This is a rather counter-intuitive result. It means that, for fixed σ and p γ , it is possible to take the system from a non-oscillating to an oscillating phase by lowering the connectivity K. 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 meanfield with quenched random graph results.

Finite-size scaling
Our phase diagrams rely heavily on the proposed order parameter q. The inset of figure 7 indicates that near σ c its behavior becomes increasingly abrupt with increasing system size N. In order to confirm that q 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 ∆ ≡ σ − σ c , FSS predicts that q ∝ L −β/ν ⊥ f ∆L 1/ν ⊥ for a lattice with linear size L [40], where the critical exponents are defined as q ∝ |∆| β and ξ ∝ |∆| −ν ⊥ in the limit N = L d → ∞ (where ξ is the correlation length). This holds for d below the upper critical dimension d c . For d ≥ d c , mean-field exponents are expected and the scaling relation has to be modified [46] with L → N 1/dc , so This modified version of the usual FSS relation was shown to also hold for infinitely coordinated networks [47] (i.e. the complete graph). For large argument, the scaling function in Equation (17) becomes f (x) ∝ x β , as usual [40]. In the subcritical regime (∆ < 0), one expects q ∼ O(N −1/2 ) [3], so f (∆N 1/dcν ⊥ ) ∝ N −1/2 N β/dcν ⊥ . For this to be true, f (x) ∝ x β− dcν ⊥ 2 when x < 0. 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 β = 1/2 [3], ν ⊥ = 1/2 and d c = 4 [40]. Similar results are obtained with complete graphs (not shown).

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 meanfield 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. [20] in a model of deterministic excitable automata (p γ = 1 in our model). The effects of setting p γ < 1, 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. [20] for strong coupling. While in their model the transition is due to very large amplitudes driving the system into rest (P t (1) = 0), in our model the system is thrown into an active albeit disordered state, which still has P t (1) = 0, 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 K.
Although we have restricted ourselves to τ = 3, preliminary results suggest that the overall scenario is preserved for larger values of τ , specially regarding the first transition at σ c . 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 N and sufficiently strong coupling, the CIC grows in amplitude and nears the absorbing state, to which the system is thrown by fluctuations [20]. 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 p γ = 1, 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. [30]. Interestingly, their model also shows an oscillating phase with reentrance, whose size decreases as memory time decreases. In our model, this corresponds to lowering p γ , 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 d c = 4 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 [48]). 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 p γ , below which no oscillations are stable. Fixing p γ < 1 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.
A u = e iθ 0 u , (A.1) with both normalized: v, u = 1 = u, u (brackets denote the standard complex inner product). Let also B( x, y) and C( x, y, z) be multilinear functions proportional to the first nonlinear terms of the Taylor expansion of F at σ = σ N S , i.e., where I is the τ × τ identity matrix and¯ u is the conjugate of u, the coefficient l 1 is finally given by [42] l 1 = 1 2 Re e −iθ 0 v, C( u, u,¯ u) + 2 v, B( u, r) + v, B(¯ u, s) .(A.7)