Excitable elements controlled by noise and network structure
We study collective dynamics of complex networks of stochastic excitable elements, active rotators. In the thermodynamic limit of infinite number of elements, we apply a mean-field theory for the network and then use a Gaussian approximation to obtain a closed set of deterministic differential equations. These equations govern the order parameters of the network. We find that a uniform decrease in the number of connections per element in a homogeneous network merely shifts the bifurcation thresholds without producing qualitative changes in the network dynamics. In contrast, heterogeneity in the number of connections leads to bifurcations in the excitable regime. In particular we show that a critical value of noise intensity for the saddle-node bifurcation decreases with growing connectivity variance. The corresponding critical values for the onset of global oscillations (Hopf bifurcation) show a non-monotone dependency on the structural heterogeneity, displaying a minimum at moderate connectivity variances.
Excitable elements are ubiquitous in physical, geophysical, chemical and biological systems. Examples range from chemical reactions, neuronal systems, cardiovascular tissues and climate dynamics to laser devices, see Ref. LiGarNeiSchi04 (); AniAstNeVaSchi07 () and references therein. An excitable system possesses a stable equilibrium, but, once excited by sufficiently strong perturbations, displays large non-monotonic excursions before returning to the state of rest.
Networks of excitable elements are omnipresent in nature, and the human brain, as a gigantic network of neurons, appears to be paradigmatic in this respect (cf. Sp10 ()). In networks of nonlinear oscillators, one typically encounters the collective phenomenon of synchronization: different oscillators adjust their rhythms due to coupling AniVadPoSaf92 (); PikRosKu03 (); BalJaPoSos10 (). Recently much effort has been invested into understanding the mechanisms which stand behind synchronization in complex networks of nonlinear oscillators and excitable systems, see Ref. ZaSaLSGNei05 (); ArDiazKuMoZh08 () and references therein.
A crucial step in characterizing limit cycle oscillators was done by Kuramoto Kur84 (): for a broad class of systems, he reduced description to dynamics of a single scalar variable, the uniformly rotating phase. Since the origins of this approach can be traced back to the biologically motivated Winfree model Win67 (), it appears to be a good candidate for providing fundamental insights into the functioning of coupled biological oscillators. Specifically, in the context of computational neuroscience various questions can be addressed by studying the Kuramoto approach with suitable modifications cumin2007 (); BreaHeiDaff10 (). In general, coupled limit-cycle oscillators can be approximated by the Kuramoto model as long as the interaction between them stays weak, so that each oscillator experiences only a small phase perturbation. Relative simplicity of the Kuramoto model turns it into a convenient tool for analytical explorations. In the frame of this model, characteristics of different transitions between desynchronized and synchronous states in large ensembles of oscillators can be obtained in closed form. This explains continuing interest in the Kuramoto model among physicists, neuroscientists and mathematicians Str00 (); AcBoViRiSp05 ().
Phase oscillators in their canonical form are not directly applicable to studies of coupled excitable elements. In particular a Kuramoto oscillator possesses translational invariance with respect to phase shifts: it rotates uniformly, and all phase values are dynamically equivalent, whereas for an excitable element a certain event (e.g. spike) is singled out. An “active rotator” – a modification of the Kuramoto model which takes into account these properties of excitable elements – was introduced by Shinomoto and Kuramoto in ShiKur86 (). Since then, the analysis was refined and certain additional aspects were investigated in KoLuLSG02 (); TesScToCo07 (); ChStr08 (); NaKo10 (); LafCoTo10 (). Most of the research considered the case of uniform global coupling in the ensemble; only in a few works, active rotators were studied on a network with complex structure, see e.g. TesZanTo08 (). Here we study dynamics of networks of randomly connected stochastic active rotators and put emphasis on how temporal fluctuations and the network structure influence the dynamics of the network as a whole.
The paper is organized as follows. In section 2 we introduce the basic model and briefly discuss its properties. In section 3 the mean-field approximation is developed: the Fokker-Planck equation for the time-dependent probability distribution of the phases is derived and transformed by means of the Fourier expansion into the infinite hierarchy of ordinary differential equations. In section 4 this hierarchy is brought to a closed form with the aid of a Gaussian approximation. Our approach is exemplified for regular and binary random networks in section 5. Numerical simulation of the binary network corroborates the conclusions drawn from the analysis of the mean-field model.
In the population of noise-driven active rotators, introduced by Shinomoto and Kuramoto ShiKur86 (), dynamics of individual phases is governed by equations
where the units are indexed by . Parameter determines the excitation threshold of each separate element, whereas characterizes the coupling strength. Here, we restrict ourselves to the case where these parameters are identical for all elements of the ensemble. The network is assumed to be undirected and unweighted, hence the elements of the symmetric adjacency matrix adopt only two values: , if the units and are coupled, and otherwise. The number of connections of a node, given by
is called the degree of this node. In absence of isolated nodes, . Complex topologies of real-world networks can be encoded into the adjacency matrix, and decoded by counting all the degrees, which gives rise to a degree distribution . Throughout this paper we discuss random networks.
The terms in (1) are statistically independent across the network nodes and are zero mean Gaussian white noise sources, where the angular brackets denote averages over different realizations of the noise and is the noise intensity.
For an isolated element in the absence of noise, , excitable behavior is apparent. For the stable equilibrium is at , and the system needs a sufficiently strong perturbation in order to make an excursion around the circle. Noise plays this role driving the system to escape from the equilibrium state . An escape event corresponds to the release of a single spike ShiKur86 (); KurSchu95 (); LiGarNeiSchi04 (). For the element shows oscillatory behavior with frequency . Notably, evolution of the phase is not uniform in this case: it is at the slowest near and at the fastest near .
3 Mean-field theory
A mean-field description can be constructed by assuming that all nodes in the network are statistically identical, so that they are fully characterized by their individual degrees SonnSchi12 (). The derivation of the mean-field description can be done by coarse-graining the network structure, namely by replacing the original network by a fully connected one with random coupling weights that mimic the actual connectivity (see Fig. 1 for illustration). Considering an uncorrelated random network Ichi04 (); SonnSchi12 (); SonnSagSchi13 () and requiring the conservation of the individual degrees, , the elements of the approximate adjacency matrix read
Inserting this into Eq. (1) yields a weighted mean-field approximation and effectively provides a one-oscillator description.
The system of coupled Langevin equations (1) can be characterized in terms of the transition probability density . Here, the vector is built from the phases of rotators at time whereas consists of their values at initial time . Components of the time-independent vector are the degrees of all nodes in the network. We assume that the initial phases are identically and independently distributed. The transition probability density obeys a Fokker-Planck equation (FPE) which describes the evolution of the ensemble from time to time .
Integration of the FPE over the phases , degrees and corresponding initial phases SonnSchi12 (), yields an evolution equation for the one-oscillator probability density , in which, however, the two-oscillator density is present as well (see Ref. CraDa99 () for the case of all-to-all coupling with more general coupling functions):
In the thermodynamic limit the ratio approaches the average degree . In the framework of the mean-field approximation, correlation of phases between any two oscillators can be discarded, In this way Eq. (4) becomes closed in , yielding an effective one-oscillator description; the price is, that it becomes a nonlinear FPE Fr05 ():
Nonlinearity enters this equation through the mean increment of the phase per unit time, , which depends on the density via the mean-field amplitude and phase ,
where denotes the degree distribution. Since all oscillators are statistically identical in this setting, from now on we neglect the indices for the individual phases and degrees. We note that Eq. (6) suggests to introduce order parameters for each degree class, namely
Let us proceed with the Fourier series expansion of the one-oscillator probability density:
with and . Through the inverse transform of (8) one can write
where ; stands for the number and for the set of nodes with degree . For every degree that appears in the network, the Fourier coefficients are governed by an infinite chain of coupled complex-valued equations
with the abbreviation . The mean degree of the network is the first moment of the degree distribution. As expected, in the case of all-to-all connectivity the results of SakShiKur88 () are recovered. Equation (10) provides a complete description of our system in the mean-field approximation, which is valid as long as . Thus, for a sparsely connected network the mean-field description is likely to fail ResOttHu05 (); SonnSchi12 (); PerErBaTirJe13 (). We also note that the above mean field approach is restricted to networks with finite second moment of the degree distribution ResOttHu05 (); Lee05 (); SonnSchi12 ().
4 Gaussian approximation
There seems to be no gain in reduction of complexity from recasting the Langevin dynamics (1) into the mean field description (10). One approach to turn (10) into a dynamical system of finite order is to truncate the Fourier series. In the case of all-to-all coupling the values of decay with increase of sufficiently fast, and the quantitative dynamics is correctly captured by the first two dozens of coefficients so that the rest of the series can be neglected SakShiKur88 (); ZaNeFeSch03 (). An alternative approach lies in finding a framework in which a closure of the infinite set of equations (10) can be achieved. The so-called Ott-Antonsen ansatz OttAnt08 () yields impressive results for deterministic ensembles of coupled Kuramoto oscillators, reducing exactly the infinite-dimensional set of equations to a low-dimensional system. Unfortunately, for the case of noisy systems the direct application of the Ott-Antonsen ansatz is impossible, and we are unaware of its appropriate modifications. Instead, we use a Gaussian approximation (GA): we assume that in every subset of oscillators with the same network degree, the distribution of the phases at every moment of time is Gaussian with the mean and the variance . The smaller the noise intensity , the more accurate is the GA KurSchu95 (). The noise should not dominate the collective dynamics, i.e. . For a recent utilization of the GA, see Ref. SonnSchi13 ().
In the thermodynamic limit, , the averaging in Eq. (9) in the GA yields
therefore all and are predetermined by and : , , etc. ZaNeFeSch03 (). Thus, the dimension of the reduced system of equations for the Fourier coefficients is twice the number of different degrees present in the network. Denoting by and by , Eqs. (10) turn into
for every degree in the network. Averaging is made over the degree distribution . Equivalently, equations for the mean and the variance are
5 Results: regular and binary random networks
The integro-differential character of Eqs. (13) makes further general analysis cumbersome. Below we discuss explicitly two relatively simple examples of regular and binary random networks.
5.1 Regular networks
For regular networks where all nodes have the same degree Eqs. (13) simplify to two coupled equations:
The only difference to the known cumulant equations for a fully connected network ZaNeFeSch03 () is a rescaling of the coupling strength: (compare panels (a) and (b) in Fig. 2). Figure 2 panel (a) shows a complete bifurcation diagram for all-to-all coupling on the parameter plane excitability and noise intensity . Two saddle-node bifurcations, denoted by and , separate parameter regions in which only a single steady state exists from the central wedge-shaped region in which three steady states are present. While at a saddle point and a node are created, destroys the saddle and the original steady state. On both bifurcation curves one of the Jacobian eigenvalues vanishes. There is a point where the second eigenvalue vanishes as well. This codimension-2 Bogdanov-Takens (BT) bifurcation is an origin of two further bifurcation lines: the Hopf bifurcation and the homoclinic bifurcation . While above , the steady state is stable, below the stable limit cycle comes into existence. The curve marks the existence of an orbit homoclinic to the saddle point. In the parameter region between the curves , and the stable steady state coexists with the attracting limit cycle. In the region between and two stable steady states exist. At the point , saddle-node and homoclinicity merge together ZaNeFeSch03 ().
Above the Hopf bifurcation line , individual units are firing incoherently, and therefore do not create a globally synchronized oscillation. So in that stationary state, even though and are constant in time, the current does not vanish ( stands for time average); it rather grows with further increase of until saturation TesScToCo07 (). Below , small-amplitude oscillations are born, where is localized between two phase values (“localized oscillations”). With decreasing noise intensity their amplitude grows. Below the line , oscillations turn into full rotations: instead of oscillating back and forth, the ensemble rotates around the whole circle . Upon the line , becomes infinite; below it is finite again. Finally, below the saddle-node bifurcation line, the system reaches again a stationary state, but this time the oscillators are mostly at rest without producing any current .
In the following, we do not discriminate between localized oscillations and full-circle rotations. We focus on how oscillating patterns are affected by the network structure. Therefore, we concentrate on finding the Hopf and the saddle-node bifurcations.
Since a reduction in the number of connections effectively decreases the coupling strength, no new qualitative phenomena are expected if all nodes share the same degree. In what follows, we use the abbreviation for the connectivity fraction.
We detect numerically the dependence on of Hopf and saddle-node bifurcations in Eqs. (14) for different excitation thresholds with the help of the software package MATCONT DhGoKuz03 (). Results are summarized in Fig. 3 showing the bifurcation value of noise intensity (called the critical noise intensity) versus for fixed values of the excitability parameter . We observe that both the Hopf and the saddle-node bifurcations are shifted to smaller noise intensities when the connectivity is decreased. Dependence of the Hopf threshold on is practically linear, in accordance with the expansion
As one would expect, the oscillatory regions shrink for higher values of threshold .
5.2 Binary random networks
A simple example of a complex heterogeneous network is a binary random network which possesses two distinct connectivity degrees. Denoting these as and (, ), the degree distribution is given by
We observe that these networks tend to be disassortative in the sense that nodes with different degrees are preferentially connected New02 (). Such degree-degree correlations could probably be factored in by developing a theory along the lines of BoPast02 (). It turns out (not shown here) that the disassortativity attains large values, if the variance among the connections is large. The maximum variance that can be obtained with (16) is . We consider only values up to , where the disassortativity is negligible.
Bifurcation analysis of this system was performed using MATCONT DhGoKuz03 ().
Fig. 4(a) illustrates a binary random network. In (b) we show the order parameters (7) obtained by simulation of the full dynamics (1). More connections appear to increase the similarity between the phases in the corresponding population. Fig. 4(c) shows phase portraits of the GA Eqs. (17) near the saddle-node bifurcation. The left panels depict the oscillatory regime. Indeed, the set of nodes with more connections is more synchronized both in the case of localized oscillations (upper left) and for full rotations (lower left). This is signalled by smaller variances . In the right panels of (c) a tiny decrease in the noise intensity brings the whole system into a stationary state, where the variables approach a fixed point.
In this context we never observed that one population of a certain degree is oscillating, while the other is not. This might be the result of entanglement of two populations with different connectivity degrees. Hence, a visualization such as in Fig. 4(a) has to be taken with care.
In order to isolate the effects of network heterogeneity, we fix the average degree, . The variance of is tuned by varying three parameters: , and . Since is not uniquely defined and different parameter sets yield the same variance, we take the path where the -values are the largest, because then the mean-field assumption is best fulfilled. Fig. 5 shows critical noise intensity at which a Hopf or saddle-node bifurcations occur in the GA system (17) versus the connectivity variance . In the oscillatory regime, , shown in Fig. 5(a), vs shows trends similar to the case of SonnSchi12 (),
Nonlinear dependence of vs starts to appear only near the threshold .
In the excitable regime [Fig. 5(b)] the oscillatory region grows with the increase of , if the variance is already large. In particular, the Hopf bifurcations shift to larger noise values, while the saddle-node bifurcations shift to smaller ones. Furthermore, for the Hopf bifurcation the dependence vs becomes non-monotone with minima at moderate values of the structural heterogeneity. These minima become more pronounced for larger values of . This effect constitutes a qualitative difference to the oscillatory regime, , and is induced by the network structure.
In particular, the minima observed for the Hopf bifurcation lines result from adding a small fraction of highly connected nodes, so-called hubs. If the nodes with larger degree are in the majority, then the critical noise is larger and the minima disappear (not shown here).
5.3 Numerical simulations of binary random networks
To test theoretical predictions made with the mean-field theory in Gaussian approximation, we performed numerical simulations of random binary networks of stochastic rotators. The goal was to check whether the global oscillations can indeed be suppressed by increasing the network heterogeneity from a small to a moderate value. We took networks of rotators and integrated Eq. (1) using the Heun method Man00 () with the time step of . Smaller time steps down to had no significant effect on the simulation results.
We used the Viger-Latapy algorithm to generate the binary random networks ViLa05 (). In particular, we used the C library igraph CsNe06 () (version 0.6.5), where the aforementioned algorithm is implemented.
Fig. 6 shows simulation results for two different networks. Both networks had the same average degree, but distinct variances. Compared with Fig. 5, the only difference in parameters of Fig. 6 is a larger value of the coupling strength which has no qualitative effect on global dynamics but simplifies the detection of the global oscillations. As can be seen from panels (a) and (b), an increase of the network heterogeneity leads to suppression of network’s oscillations, in agreement with the theory. Panel (c) shows snippets of the phase distribution with the parameters from panel (a). As in Fig. 4(c)[upper left], we observe localized oscillations. One can recognize the distinct single-hump shape that is necessary for the Gaussian approximation (section 4).
We have studied a well-known model for excitable dynamics, the so-called active rotators, in the context of heterogeneous networks. Besides investigating the effects of noise, we put emphasis on the effects that emerge due to different coupling structures. Considering the infinite system-size limit, the first step in our analysis has been to reduce the dimensionality of the system. Coarse-graining of the network has lead us to the mean-field description. Assuming further that the phases in each set of oscillators with the same degree obey a Gaussian distribution with time-dependent cumulants, we arrived at a lower-dimensional set of integro-differential equations. The order of this system equals the doubled number of different degrees occuring in the network. The approach has been used to perform numerical bifurcation analysis on regular and binary random networks. The calculated state diagrams display the parameter regions where the neuronal network shows synchronous or asynchronous spiking behavior. Transitions between these regions are typically related to passages through the Hopf and saddle-node bifurcations. When the noise intensity is increased from zero, the saddle-node bifurcations mark the onset of global firing in the ensemble, whereas the Hopf bifurcation leads to the disappearance of the collective mode of synchronous oscillations. We have been able to show that changing the network structure alone, without variations in the noise intensity and/or the excitation threshold, is sufficient for a switch between the different patterns of spiking behaviors. In fact, adaptations in the coupling structure play a crucial role for the information processing in populations of neurons Sp10 ().
Our results suggest that deleting uniformly connections of all nodes in the network does not lead to new characteristics, since the coupling strength is only rescaled effectively. However, if nodes differ in their number of connections, an interesting phenomenon is obtained in the excitable regime. We have found out that the critical noise intensity, at which the Hopf bifurcation occurs and small-scale oscillations are born, can possess a minimum for a moderate network heterogeneity. It was induced by the presence of a small fraction of highly connected nodes. This might have meaningful implications for the functioning of neuronal networks, since such hub nodes naturally exist Sp10 (); VarChPanHaChk11 (). It is well-known that neuronal networks are highly heterogeneous. Our analysis further suggests that this property implies existence of large parameter regions to get global firing. If, however, a specific task requests a switching between different spiking patterns, reducing the heterogeneity in the connectivities becomes beneficial, since only small noise variations are needed. The optimum for such a task would then be given by a moderate coupling variance with a small fraction of hub nodes.
Acknowledgements.This work was supported by the GRK1589/1 and project A3 of the Bernstein Center for Computational Neuroscience Berlin. Research of M.Z. was supported by the project D21 of the DFG Research Center MATHEON.
- (1) B. Lindner, J. García-Ojalvo, A. Neiman, and L. Schimansky-Geier. Phys. Rep., 392:321, 2004.
- (2) V. S. Anishchenko, V. Astakhov, A. Neiman, T. Vadisova, and L. Schimansky-Geier. Nonlinear Dynamics of Chaotic and Stochastic Systems. Springer-Verlag, Berlin, 2007.
- (3) O. Sporns. Networks of the Brain. MIT Press, 2010.
- (4) V. S. Anishchenko, T. E. Vadivasova, D. E. Postnov, and M. A. Safonova. Int. J. Bifurcation Chaos, 02:633, 1992.
- (5) A. Pikovsky, M. Rosenblum, and J. Kurths. Synchronization: A universal concept in nonlinear sciences. Cambridge Univ. Press, U. K., 2003.
- (6) A. Balanov, N. Janson, D. Postnov, and O. Sosnovtseva. Synchronization: From Simple to Complex. Springer-Verlag, Berlin, 2010.
- (7) M. A. Zaks, X. R. Sailer, L. Schimansky-Geier, and A. B. Neiman. Chaos, 15:026117, 2005.
- (8) A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou. Phys. Rep., 469:93–153, 2008.
- (9) Y. Kuramoto. Chemical Oscillations, Waves, and Turbulence. Springer-Verlag, Berlin, 1984.
- (10) A. T. Winfree. J. Theor. Biol., 16:15–42, 1967.
- (11) D. Cumin and C. P. Unsworth. Physica D, 226(2):181–196, 2007.
- (12) M. Breakspear, S. Heitmann, and A. Daffertshofer. Front. Hum. Neurosci., 4:190, 2010.
- (13) S. H. Strogatz. Physica D, 143:1–20, 2000.
- (14) J. A. Acebrón, L. L. Bonilla, C. J. Pérez Vicente, F. Ritort, and R. Spigler. Rev. Mod. Phys., 77:137–185, 2005.
- (15) S. Shinomoto and Y. Kuramoto. Prog. Theor. Phys., 75(5):1105, 1986.
- (16) M. Kostur, J. Luczka, and L. Schimansky-Geier. Phys. Rev. E, 65:051115, 2002.
- (17) C. J. Tessone, A. Scirè, R. Toral, and P. Colet. Phys. Rev. E, 75:016203, 2007.
- (18) L. M. Childs and S. H. Strogatz. Chaos, 18:043128, 2008.
- (19) K. H. Nagai and H. Kori. Phys. Rev. E, 81:065202(R), 2010.
- (20) L. F. Lafuerza, P. Colet, and R. Toral. Phys. Rev. Lett., 105:084101, 2010.
- (21) C. J. Tessone, D. H. Zanette, and R. Toral. Eur. Phys. J. B, 62:319, 2008.
- (22) C. Kurrer and K. Schulten. Phys. Rev. E, 51(6):6213, 1995.
- (23) B. Sonnenschein and L. Schimansky-Geier. Phys. Rev. E, 85:051116, 2012.
- (24) T. Ichinomiya. Phys. Rev. E, 70:026116, 2004.
- (25) B. Sonnenschein, F. Sagués, and L. Schimansky-Geier. Eur. Phys. J. B, 86:12, 2013.
- (26) J. D. Crawford and K. T. R. Davies. Physica D, 125:1–46, 1999.
- (27) T. D. Frank. Nonlinear Fokker-Planck Equations. Springer-Verlag, Berlin Heidelberg, 2005.
- (28) H. Sakaguchi, S. Shinomoto, and Y. Kuramoto. Prog. Theor. Phys., 79(3):600, 1988.
- (29) J. G. Restrepo, E. Ott, and B. R. Hunt. Phys. Rev. E, 71:036151, 2005.
- (30) T. Pereira, D. Eroglu, G. B. Bagci, U. Tirnakli, and H. J. Jensen. Phys. Rev. Lett., 110:234103, 2013.
- (31) D.-S. Lee. Phys. Rev. E, 72(026208):1–6, 2005.
- (32) M. A. Zaks, A. B. Neiman, S. Feistel, and L. Schimansky-Geier. Phys. Rev. E, 68:066206, 2003.
- (33) E. Ott and T. M. Antonsen. Chaos, 18:037113, 2008.
- (34) B. Sonnenschein and L. Schimansky-Geier. Phys. Rev. E, 88:052111, 2013.
- (35) A. Dhooge, W. Govaerts, and Y. A. Kuznetsov. ACM TOMS, 29:141–164, 2003.
- (36) M. E. J. Newman. Phys. Rev. Lett., 89(20):208701, 2002.
- (37) M. Boguná and R. Pastor-Satorras. Phys. Rev. E, 66:047104, 2002.
- (38) R. Mannella. A Gentle Introduction to the Integration of Stochastic Differential Equations. In J. A. Freund and T. Pöschel, editors, Stochastic Processes in Physics, Chemistry, and Biology, volume 557, pages 353–364. Springer-Verlag, Berlin Heidelberg, 2000.
- (39) F. Viger and M. Latapy. Efficient and Simple Generation of Random Simple Connected Graphs with Prescribed Degree Sequence. In L. Wang, editor, Comput. Combin., volume 3595, pages 440–449. Springer Berlin Heidelberg, 2005.
- (40) G. Csárdi and T. Nepusz. InterJournal Complex Systems, (1695), 2006.
- (41) L. R. Varshney, B. L. Chen, E. Paniagua, D. H. Hall, and D. B. Chklovskii. PLoS Comput. Biol., 7(2):e1001066, 2011.