- Email: [email protected]
Type I vs. type II excitable systems with delayed coupling
Chaos, Solitons and Fractals 23 (2005) 1221–1233 www.elsevier.com/locate/chaos
Type I vs. type II excitable systems with delayed coupling Nikola Buric´
a,*
, Ines Grozdanovic´ b, Nebojsˇa Vasovic´
b
a
b
Department of Physics and Mathematics, Faculty of Pharmacy, University of Beograd, Vojvode Stepe 450, Beograd, Serbia and Montenegro Department of Applied Mathematics, Faculty of Mining and Geology, P.O. Box 162, Beograd, Serbia and Montenegro Accepted 8 June 2004
Abstract Excitable and oscillatory dynamics of delayed locally coupled type I and type II excitable systems is analyzed. Diffusive and sigmoid coupling have been considered. It is shown that the stability and the patterns of exactly synchronous oscillations depend on the type of excitability and the type of coupling. However, within the same class, characterized by the excitability type and the coupling the dynamics qualitatively depends only on whether the number of units is even or odd. 2004 Elsevier Ltd. All rights reserved.
1. Introduction Excitability is a common property of many physical and biological systems. Although there is no precise definition [1] the intuitive meaning is clear: A small perturbation from the single stable stationary state can result in a large and long lasting excursion away from the stationary state before the system is returned back asymptotically to equilibrium. Further more, as an external parameter is changed, the global attractor in the form of the stationary point bifurcates into a stable periodic orbit, and the excitability is replaced by the oscillatory dynamics. The are basically two qualitatively different types of excitability, fancifully called type I and type II, which are distinguished by different bifurcations leading to the oscillatory dynamics [2]. This is manifested in different properties of the frequency of the emergent oscillations. Our goal in this paper is to study the differences in the dynamics of relatively small chains of coupled type I or type II excitable systems. As the example of the type I excitability we shall use the Terman–Wang model [3] for the corresponding values of the parameters, and as the example of the type II excitable system we shall use the FitzHugh–Nagumo model [4,5]. Both systems originate as models or neuronal excitability, i.e. of the dynamics of the neuron transmembrane potential, but this interpretation of the models will not be of primary importance for us. We should stress that we shall analyze the systems such that each of the units in a chain is a system operating in the excitable, and not in the oscillatory, regime. The oscillations might appear only due to the coupling between the excitable units. Transfer of excitations between the neighboring units is not instantaneous, so that realistic models of the coupling should depend on the states of the neighboring units at different times. Thus we shall suppose that the function which
*
Corresponding author. E-mail address: [email protected] (N. Buric´).
0960-0779/$ - see front matter 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2004.06.033
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
1222
describes the influence of the ith unit on the jth unit at the time t depends on the state of the ith unit at some earlier time t s. In fact we shall consider two types of coupling: the so-called diffusive, described by f(xj, xi) = xj(t) xi(t s) and the sigmoid type given by f ðxj ; xi Þ ¼ tan1 xi ðt sÞ. The effects of time-lags on the dynamics of coupled excitable systems are much less studied than for the systems of delayed coupled regular or chaotic oscillators. Let us illustrate typically studied systems by quoting only some representative references. Delayed coupled relaxation oscillators have been studied for example in [6–8]. In these studies, each unit is a relaxation oscillator, and the primary objective of the analysis are the periodic orbits that appear in the delayed coupled system. Singular approximation, or an approximate or numerical construction of the Poincare map, are used to analyze various types of synchronous or asynchronous oscillations. In [9,10], a network of N = 2 and N > 2 oscillators described by the equations of the normal form of the Hopf bifurcations with delayed diffusive coupling is studied. It is shown that the time-delay can lead to the stabilization of the trivial stationary point, which is interpreted as the amplitude death. Less directly related to our work is the analyzes of the influence of the time delay in the systems of coupled phase oscillators. The influence of time-delay in coupled phase oscillators was studied for example in [11,12], where it is shown that the time-delay cannot introduce significant changes into the dynamics of a class of such systems, unless the time-lag is of the order of several oscillation periods [12]. Independently of neuronal models, collective behavior of the phase coupled (phase) oscillators with time-delayed coupling, have been studied using the dynamical [13–15], or statistical [16] methods. Finally, the influence of time delay has been studied in the Cohen– Grasberger–Hopfield type of neural networks, as early as in 1967 [17]. More recent references are for example [18,19], and for small networks [20,21] (see also [22] and the references there in). Delayed coupled type II excitable, as opposed to oscillatory, systems have been recently studied in reference [23] (where a collection of references to the typical results for coupled oscillators is given). This was extended in [24] to chains of FitzHugh–Nagumo systems with small in-homogeneities. In these papers only the type II system, given by the FitzHugh–Nagumo equations, and the sigmoid delayed coupling is considered. In this paper we shall compare the characteristic features of delayed coupled type I and type II systems for different typical forms of the coupling, and consider in more details the properties of the patterns of synchronous oscillations that occur in these systems. The paper is organized as follows. In Section 2 the two basic models of type I and type II systems and the two types of coupling will be introduced. In Section 3 the local bifurcations of the relevant stationary solutions will be analyzed. Stable oscillatory dynamics will be studied in Section 4. Section 5 contains the summary of our main conclusions.
2. Two basic models In this section we introduce the class of models that will be analyzed. We shall consider one-dimensional chains of N excitable units with periodic boundary conditions, where only the nearest neighbors are coupled. Such a system is given by the equations of the following general form: x_ i ¼ X ðxi ; y i Þ þ Cðxi ; xiþ1 ðt sÞ; xi1 ðt sÞÞ; y_ i ¼ Y ðxi ; y i Þ; i ¼ 1; 2; . . . ; N ;
ð1Þ
with the following boundary conditions on Nth and 0th unit: ðx0 ; y 0 Þ ¼ ðxN ; y N Þ;
ðxN þ1 ; y N þ1 Þ ¼ ðx1 ; y 1 Þ:
ð2Þ
The field X(x, y), Y(x, y) describes the single isolated unit, and the function Cðxi ðtÞ; xii ðt sÞ; xiþ1 ðt sÞÞ Cðxi ðtÞ; xsi1 ðtÞ; xsiþ1 ðtÞÞ describes the time-delayed coupling between the nearest neighbors. As the examples of type I and type II excitability we shall take the Terman–Wang and FitzHugh–Nagumo models, respectively. In the first case, the single unit is described by X ðx; yÞ ¼ X TW ¼ x3 þ 3x y þ 2 þ I; Y ðx; yÞ ¼ Y TW ¼ ½að1 þ tanhðx=bÞÞ y;
ð3Þ
where we fix the parameters to the following typical values a = 3, b = 0.2, I = 0.4, = 1, so that the system has a stable stationary solution (x0, y0) (1.34577, 0), which would go through the saddle-node on the limit cycle bifurcation upon increasing the absolute value of the external parameter I. The resulting stable limit cycle has finite radius and the frequency is vanishingly small. This type of bifurcation characterizes type I excitability. For type II excitability we shall take the FitzHugh–Nagumo model in the form [25]
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
1223
Fig. 1. The nulclines and a typical excitable orbit for the type I system (). The parameters are a = 3, b = 0.2, I = 0.4 and the stationary state is at (x0, y0) (1.34577, 0).
X ðx; yÞ ¼ X FN ¼ x3 þ ða þ 1Þx2 ax y þ I; Y ðx; yÞ ¼ Y FN ¼ bx cy:
ð4Þ
where a, b and c are positive internal parameters, and I, which will be taken I = 0, is the external parameter. The stationary solution (0, 0) is stable for 4b < c(a 1)2, I = 0. The stationary solution loses stability through the Hopf bifurcation upon increasing I and a stable limit cycle with vanishingly small radius but finite frequency is born. This type of bifurcation characterizes type II excitability. Notice that x_ nullcline in both models have the same quibic shape, and the main difference between the models is in y_ nullcline, which is linear for the FitzHugh–Nagumo and sigmoid for the Terman–Wang. The nullclines and typical excitable orbits are presented in Figs. 1 and 2, for the Terman–Wang and the FitzHugh–Nagumo systems, respectively. Two types of coupling will be analyzed. The diffusive coupling, Cðxi ðtÞ; xsi1 ðtÞ; xsiþ1 ðtÞÞ ¼ cð2xi þ xsi1 þ xsiþ1 Þ
ð5Þ
and the sigmoid coupling in the following form: Cðxi ðtÞ; xsi1 ðtÞ; xsiþ1 ðtÞÞ ¼ ctan1 ðxsi1 x0 Þ þ ctan1 ðxsiþ1 x0 Þ;
ð6Þ
where (x0, y0) is the stable stationary solution. Thus, the coupling between the units in equilibrium is zero. As expected, the model of the coupling plays important role in the dynamics of the whole chain, but, nevertheless, some properties of the synchronization patterns are independent of the type of the coupling.
3. Bifurcations of stationary states In this section we present results of the analyses of local bifurcations due to time-lag of the stationary solution for a system of only two delayed-coupled excitable units. Only the stationary solution which corresponds to the excitable dynamics will be considered. Coordinates of the stationary solution do not depend on the coupling constant c or the time-lag s, but the stability could be effected by varying either c or s. We shall be interested in the dependence
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
1224
Fig. 2. The nulclines and a typical excitable orbit for the type II system (). The parameters are a = 0.25, b = c = 0.02, and the stationary state is at (x0, y0) = (0, 0).
of the bifurcation values of the time-lag s on the coupling constant c, for all intrinsic parameters of each unit fixed to some typical values, such that the considered stationary solution is stable when c = 0 and s = 0. The bifurcations occur for the values of the parameters and s such that there is a non-hyperbolic root, k, Re k = 0 of the characteristic equation. As we shall see the relevant characteristic equation has the same general form for all considered models and couplings, the differences between the models occurring in the explicit relations between the parameters and the coordinates of the considered stationary solutions. Linearization of the any of the considered models, around an arbitrary stationary solution x10 = x20 ” x0 and y10 = y20 ” y0 and substitution: xi(t) = Ai exp kt; yi = Bi exp kt results in the following characteristic equation: DðkÞ ½ðA þ kÞð þ kÞ þ B2 ½cð þ kÞ expðksÞ2 ¼ 0;
ð7Þ
which can be factored as D(k) = D1(k)D2(k) where D1;2 ðkÞ ¼ ½k A þ B c expðksÞð þ kÞ þ B ¼ 0:
ð8Þ
In the Terman–Wang case we have A ¼ 3 3x20 ;
B¼
a bcosh2 ðx0 =bÞ
ð9Þ
for the sigmoid coupling, and for the diffusive coupling A is replaced by A c. In the FitzHugh–Nagumo case we have A a;
B b;
¼c
ð10Þ
for the sigmoid coupling, and for the diffusive coupling A is replaced by A c, analogously as for the type I case. The difference between the systems and the couplings is further manifested in the characteristic equation when the explicit expressions for the respective stationary solutions (x0, y0) are substituted. Let us consider first the zero solutions of (7). This type of bifurcations occurs also when s = 0. Substitution k = 0 implies
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
A B c ¼ 0:
1225
ð11Þ
Next we consider pure imaginary solutions of D(k) = 0. If some additional non-degeneracy conditions are satisfied this corresponds to the Hopf bifurcation. Substitution k = ix where x is assumed real and positive, into D1 = 0 or D2 = 0, and separation of the imaginary and real parts of the resulting equations gives a biquadratic equation for x, x4 þ Mx2 þ N ¼ 0;
ð12Þ
where M ¼ A2 þ 2 2B c2 ;
N ¼ A2 2 þ B2 2AB c2 2
with the formal solutions: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2 ¼ ðM M 2 4NÞ=2;
ð13Þ
ð14Þ
which of the solutions x± are real depends on the model and its stationary solution. For the stationary solutions which represent the quiescent states of the units, we shall see that, only the pair of FitzHugh–Nagumo systems with the sigmoid or the diffusive coupling do have both type of solutions, i.e. both x2þ and x2 are real and positive. In the case of the type I systems with either types of coupling, the only real solution, for the corresponding unexcited stationary state, is x+. The critical values of the time-lag sc as a function of the parameters is obtained by inverting: tanðsc Þ ¼
xð2 þ x2 BÞ ; A2 Ax2 B
ð15Þ
where the upper or lower signs correspond to D1/2 = 0, and substituting the real and positive solutions for x of (12) in (15). The nature of the Hopf bifurcation for s = sc is determined by the dependence of the real part of k on s as it goes through sc. This is given by the sign of d Re k=ds for s = sc. Using the factorized characteristic equation, one obtains d Re k A2 þ 2 2B þ 2x2 c2 ¼ : ð16Þ ds s¼sc c2 ðx2 þ 2 Þ The previous expression is positive for x = x+, and the Hopf bifurcation increases the number of unstable direction as s is increased through sc. In the case there is real and positive x the expression (16) is negative for x = x, and the Hopf bifurcation decreases the number of unstable directions as s is increased through sc. As we have pointed out, the latter case occurs only for the FitzHugh–Nagumo systems. In what follows we denote by scj; , j = 0, 1, 2, . . . , different branches of the bifurcation curve given by (14) and (15). Let us now apply the formulas (11), (14) and (15) to obtained the bifurcations curves in the (c, s) plane for particular models and the relevant stationary solutions. Thus, we fix the internal parameters of each unit to such values that the stationary solution which correspond to the unexcited state is stable. The values of the coordinates of the stationary solution (x0,y0) do not depend on the coupling constant c. Internal parameters and the corresponding (x0, y0) are substituted into the formulas for x and sc. The result is the set of bifurcation curves scj; for the considered stationary solution. 3.1. Type I systems In the case of the Terman–Wang model, we take for both units = 1, a = 3, I = 0.4, b = 0.2. The point (x1, 0, y1, 0, x2, 0, y2, 0) = (1.34577, 0, 1.34577, 0) approximates the stationary solution up to the fifth decimal place. The bifurcation curves for the sigmoid and the diffusive coupling are presented in Fig. 3a and b. We point out that the solution x+ is positive only for such values of the coupling constant c for which the stationary solution is unstable when s = 0. There is a set of destabilizing Hopf bifurcation curves scj;þ in the (c, s) plane. The solution x is always negative, so there are no stabilizing Hopf bifurcations. Thus, the local effect of the time-delay on the stationary solution is only an increase of the number of unstable directions of the considered stationary solution. On the other hand, stable periodic orbits could appear due to bifurcations of other stationary solutions, which are stable for s = 0 and the considered c, or due to some global bifurcation. Thus, the existence of stable oscillatory solutions is not indicated by the local analysis of the considered stationary solution, and will be analyzed by numerical calculations, as is presented in the next section.
1226
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
Fig. 3. Hopf bifurcation curves scj;þ ðcÞ for the type I system and (a) sigmoid (b) diffusive coupling. The fixed values of the local parameters are as in Fig. 1.
3.2. Type II systems The Hopf bifurcation curves for two diffusively or sigmoid coupled FitzHugh–Nagumo, i.e. type II systems, are presented in Fig. 4. The values of the internal parameters are for both units a = 0.25, b = c = 0.02. The stationary solution is (x1, 0, y1, 0, x2, 0, y2, 0) = (0, 0, 0, 0) and, due to the polynomial form of the equations, it is found exactly. The stationary solution is stable for these values of the internal parameters and jcj 6 c0. The case of sigmoid coupled FitzHugh–Nagumo units was thoroughly discussed in [23,24]. Important difference, as compared with type I systems, is that even the instantaneous coupling of excitable type II units leads to a supercritical Hopf bifurcation at c0 = (a + c)/d, c < 1, which results in a stable periodic orbit. For the sigmoid coupling, the dynamics of the two units on the resulting periodic orbits is exactly synchronized, i.e. x1(t) = x2(t); y1(t) = y2(t). Further important characteristic feature of this case is the existence of real x solutions and the corresponding stabilizing Hopf bifurcations on scj; , j = 1, 2, . . . As was discussed in the cited references, this together with a global fold limit cycle bifurcation leads to the phenomenon of the oscillator death in a certain domain of the (c, s) plane. The diffusive coupling has similar effect on the stability of the stationary solution. Again, there is a supercritical Hopf bifurcation due to the instantaneous coupling (s = 0), at some c = c0 = (ac + c) < 0, and the stationary solution becomes unstable. However, the stable periodic orbit is born in the plane x1 = x2;y1 = y2 and corresponds to the out of phase oscillations. Non-zero time-lag can lead to either stabilizing or destabilizing Hopf bifurcations. In fact, both solu-
Fig. 4. Hopf bifurcation curves scj;þ ðcÞ for the type II system and (a) sigmoid (b) diffusive coupling. The fixed values of the local parameters as in Fig. 2.
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
1227
tions x+ and x can be real. Again, the oscillator death occurs in a certain domain of the (c, s) plane. Let us stress again that the phenomenon of oscillator death does not occur for coupled type I systems.
4. Oscillatory dynamics Our goal in this section is to describe the main types of oscillatory solutions and to study the stability of synchronization between different units. Finally, we shall compare the conclusions about type I systems with the corresponding behavior of type II systems. In order to describe the stable periodic solutions of N coupled units we have calculated numerically the mean frequencies of the projections of the particular stable periodic orbit on (xi, yi) plane. Stability of the exact synchronization is studied by calculation of the largest transversal Lyapunov exponent. Let us briefly recapitulate the definitions of the mean frequency and the transversal Lyapunov exponent for DDE. The mean frequency of the ith unit is introduced via the corresponding phase which is commonly defined by [26,27] ui ðtÞ ¼ 2p
t tki þ 2pk; tki
tkþ1 i
ð17Þ
where tki is the time of the kth crossing of xi(t) of some convenient threshold. If the intersection times are eventually equally distributed then the curve ui(t) is asymptotic to a straight line, i.e. limt!1 ui ðtÞ=xi t ¼ 1. The constant xi is the mean frequency, with the corresponding period Ti = 2p/xi, of the ith unit. The concept of the mean frequency enables us to talk about coherent dynamics even when the oscillations of each unit are not harmonic. The dynamics of two units i and j is coherent if the mean frequencies are equal, i.e. xi = xj Two coherent units i and j oscillate out-of-phase if tki ¼ tkj þ px1 . If tki ¼ tkj the two units oscillate in-phase. In order to study the stability of the exact synchronization between ith and jth unit we numerically calculate the dynamics: Di;j ðtÞ ¼ xi ðtÞ xj ðtÞ:
ð18Þ
A stable trivial stationary solution Di, j = 0 corresponds to the stable exactly synchronous solution xi(t) = xj(t). It is useful (see for example [28,29] to consider the dynamics of Di, j as a dynamical system on the phase space given by continR0 uous functions D defined on the interval [s, 0], with the norm kDk2 ¼ s D2 ðhÞ dh. Solutions for different initial functions in C[s, 0] generate a semi-flow on this phase space given by Dt(h) = D(t + h), t 2 R+, h 2 [s, 0]. The largest transversal Lyapunov exponent for the (i, j) exact synchronization is then defined in analogy with the finite-dimensional systems by 1 kDi;j;t k Li;j ¼ lim ln ; ð19Þ t!1 t kDi;j;0 k where the norm iDi,j;0i is small. (i, j) synchronization is locally attracting for some (c, s) if Li, j is negative for at least some sufficiently small Di,j,0. 4.1. Type I systems Chains of type I systems with delayed diffusive or sigmoid coupling have quite different global behavior. We shall first discuss the case of the diffusive coupling. The main results of our computations in this case are presented in Figs. 5 and 6. The parameters a, b and I are fixed to some typical values (for example, the same as in the previous section). Many initial functions are evolved, and the transients are discarded. There are no oscillatory solutions for any c and s = 0 or sufficiently small. However, stable oscillatory behavior appears for the time-lag greater then some scr. This scr and the oscillatory solutions are in no way related to the Hopf bifurcations dynamics of the considered stationary solution which was studied in the previous section. In fact, oscillatory solutions are possible only for c such that the stationary solution is unstable and scr is larger then s0,+. Typical dependence of the mean frequency x1, of the first component x1 on the stable limit cycles, on the time-lag, for N = 5 and N = 6, is illustrated in Figs. 5a and b. In fact, we found that the mean frequencies of all components are equal, so that the oscillations are coherent for any coupling c and time-lag s. Thus, Figs. 5 represent the frequency of any of the components. Furthermore, the typical dependence of the common frequency on the time-lag is approximately independent of N. The general conclusion is that for c < 0 the mean frequency decreases as s is increased. For c > 0 and s larger than some critical value scr the frequency first strongly depends on s in a small interval beyond scr.
1228
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
Fig. 5. Mean frequency xi as a function of s for diffusively coupled type I systems: (a) N = 5 and (b) N = 6. The fixed values of c are (1) c = 1 (boxes), (2) c = 4 (crosses), (3) c = 4 (triangles). The local parameters are as in Fig. 1.
Fig. 6. Coherent oscillations for diffusively coupled N = 6 type I systems: (a) out of phase (xi, xi + 1) and exactly synchronous (xi, xi + 2) for c = 1, s = 3; (b) convergence of the computations of the transversal Lyapunov exponents for c = 4, s = 20; (c) oscillations for typical parameters; (d) transient dynamics for large s = 50.
Further away from scr, the frequency oscillates with s in a quite narrow band. For small variable c, and fixed s, the frequency increases with c, and then weakly depends on the larger c. Let us stress that the same behavior is observed for any N as large as N = 50. The patterns of stable synchronization depend on N being odd or even. For any odd N the nearest neighbors always oscillate in exact synchrony, which is the simplest spatially uniform pattern. This is also the only patters that we have observed for any even N and c < 0. However, for even N, and c > 0, next to the nearest neighbors stable exactly synchronous oscillations are possible, i and i + 1 unit could oscillate out of phase, as it is illustrated in Fig. 6a. The out of
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
1229
phase dynamics between the nearest neighbors was observed only in a narrow interval of parameters c and s scr(c). In general exact synchronization of the nearest neighbors is much slower than the synchronization between the next to the nearest units. This is clearly shown by calculations of the transverse Lyapunov exponent Li;j ¼ limt!1 Li;j ðtÞ, illustrated in Fig. 6b. Let us point out that, for larger time-lags, say s = 50 like in Fig. 6d, the approach to the periodic oscillations on the attractor, could involve quite long transient dynamics, during which a sequence of almost synchronous bursts of oscillations alternate periods of stationary behavior. One could contemplate a source of a bursting behavior, which is executed by appropriate variations of the time-lag in the transmission of interaction through the junctions (synapses) between the units. The behavior of chains with type I units and delayed sigmoid coupling like in (6) is trivial. The symmetric stationary solution, that corresponds to the equilibrium, remains stable for any s, as long as the coupling c is such that the state is stable for s = 0. Larger values of c destabilize this state, new stable stationary states are introduced, but no oscillations occur. 4.2. Type II systems Instantaneous diffusive coupling can induce oscillations in a chain of type II excitable systems via Hopf bifurcation for the coupling parameter c = c0 < 0. There are no oscillations for c > 0. The resulting oscillations, for c < c0 < 0, of each unit are coherent. The transition from regular monochromatic oscillations of each unit, for c c0 to relaxation oscillations occurs in a very small interval of c < c0. Nevertheless the mean frequency of different units, as defined in (17), is the same for any c. Let us point out that if the parameter c is such that the units are excitable, i.e. the only attractor is the stable stationary solution, and c is far away from c0, then the time-lag has no effect on the qualitative properties of the dynamics. The stationary solution remains stable for any s. There are two distinct domains of the coupling c for which the dynamics is oscillatory. Firstly, there is a small interval near c0 for which the instantaneously coupled units are excitable but the time-lag does introduce stable oscillations, as it is indicated by the bifurcations curves in Fig. 4. The other situations occurs for jcj > jc0j, i.e. when the instantaneously coupled units oscillate, and when further Hopf bifurcations are possible due to non-zero time-lag. The results of our numerical computations are presented in Figs. 7–9. Fig. 7a and b, illustrate the dependence of the mean frequency on the time-lag for the values of c very close to c0, and for some typical value further away from c0. The Hopf bifurcation of the stationary solution at c0 (for s = 0) is responsible for the appearance of the oscillatory dynamics, so that the frequency near the bifurcation has a finite non-zero value. The dependence of xi(s) for fixed c is piecevise monotonically decreasing. The bifurcation values of s, corresponding to the jumps in xi(s), also correspond to changes of the patterns of synchronization between the units, i.e. to the threshold of the nearest neighbor exact synchronization. Possible patterns are illustrated in Figs. 8 and 9. The jumps do not occur if the fixed value of c is very close to c0. Again, the chains with odd or even N behave differently. For odd N and jcj > jc0jclose to, or moderately away from c0, and for small s, there is no stable exact synchronization between any two units. Larger s leads to stable exact synchronization. However, the only pattern of the stable exact synchronization is the most symmetric one, i.e. between all of the units. The synchronization threshold of the time-lag s decreases with increasing the absolute value of c < c0. For example, for c = 0.1 the threshold is s = 22 and for c = 0.3 the threshold is s = 9. Sufficiently large c induces exactly
Fig. 7. Mean frequency xi as a function of s for diffusively coupled type II systems: (a) N = 5 and (b) N = 6. The fixed values of c are (1) c = 0.1 (boxes), (2) c = 0.3 (crosses), (3) c = 0.056 (triangles). The local parameters are as in Fig. 2.
1230
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
Fig. 8. Illustrates transition from asynchronous coherent to exactly synchronous (i, i + 1) pattern as s is increased, for diffusively coupled N = 5 type II systems (). Frame 8f shows the transversal Lyapunov exponent for s = 20 (at the synchronization threshold) and s = 22 (slightly above the threshold) The local parameters are as in Fig. 1, c = 0.1.
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
1231
Fig. 9. Illustrates bifurcations of the synchronization patterns for diffusively coupled N = 6 type II systems (): (a) increasing c for fixed s = 20 leads to (xi, xi + 1) exact synchronization; (b) illustrates (xi, xi + 2) exact synchronization for s = 10 and c = 0.1.
synchronous oscillations with arbitrary small s. The synchronization threshold is clearly indicated by the transversal Lyapunov exponent. Chains with even N display two different patterns of stable exact synchronization. For the coupling c smaller than c0 but very close, the next to the nearest neighbors exactly synchronous oscillations are stable. The nearest neighbors oscillate out of phase. This pattern persists for any time-lag s. On the other hand, if the absolute value of c < c0 is larger, (i, i + 2) pattern is substituted by the simplest pattern with all units exactly synchronized. It is interesting to observe that for moderate values of s there could coexist locally stable attractors that correspond to different patterns of exact synchronization. For example, for N = 6, c = 0.1 and s = 15 we have found that the (i, i + 2) pattern is locally stable, but also the pattern with x1 = x4 and x6 = x3 = x2 has a non-zero domain of stability. Larger time-lag again leads to globally stable (i, i + 1) pattern of synchronization. Finally, we would like to make just a few comments on the oscillations and the patterns of synchronization in chains with sigmoid delayed coupling. Oscillations in chains of type II excitable units with such coupling have been analyzed before in [23,24]. Here we would like just to point out that the main differences compared with the diffusive coupling is that for zero time-lag sufficiently strong sigmoid coupling leads to stable exactly synchronous oscillations of all units. Thus we have to consider two different domains of coupling in analyzing the influence of the time-lag. If the instantaneously coupled units have excitable, and not oscillatory dynamics, large time delay can lead to a global fold limit cycle bifurcation, with a bi-stability as the result. The stationary state and a limit cycle are the two attractors. The oscillations between the neighbors are out of phase but the next to the nearest units oscillate in exact synchrony. Sufficiently strong instantaneous sigmoid coupling leads to exactly synchronous oscillations of all type II units, but already relatively small time-lags introduce the next to the nearest neighbor pattern for even N. As we have seen, the diffusive coupling leads to the opposite situation. Furthermore, the oscillations introduced by instantaneous coupling of type II systems, can be turned of by the time-delay. There is no such phenomenon in the chains of type I systems with either of the considered types of coupling.
5. Summary and discussion We have studied dynamics of chains of bidirectionally and locally coupled excitable systems. Sigmoid and diffusive coupling, with explicitly included time-delay, have been considered. Our main goal was to compare the behavior of type I vs. type II excitable systems with respect to the following two questions. The first one was the change of stability, with respect to coupling strength and time-delay, of the state that represent the attractor of the excitable dynamics. The second question was concerned with possible spatial patterns of exact synchronization between the oscillations of each unit that originate only through coupling. Let us stress that the isolated units are always only excitable and any oscillatory behavior is introduced by coupling. Concerning the first question, we have concluded that even the instantaneous coupling of either type can destabilize the stationary solution of the excitable behavior. However, for type II systems the destabilization occurs via Hopf bifurcation and a small stable limit cycle appears in a neighborhood of the considered stationary solution. On the other hand, the destabilization of the type I stationary state by the instantaneous interaction just introduces new stable stationary solution. In correspondence with this difference between the two types of excitable systems, the time-delay has
1232
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
different effect on the dynamics near the stationary solution. For type I systems the time-delay can only increase the number of unstable directions via supercritical Hopf bifurcation. For type II systems the time-delay can also decrease the number of unstable directions via inverse subcritical Hopf bifurcation. In particular, the time-delay can induce the oscillators death only in the chains of type II systems. Concerning the second question, we first point out that oscillatory solutions in chains of type I excitable systems are observed only for non-zero time-lag and diffusive coupling. In this case and for odd N, the only type of spatial pattern of exact synchronization is the most symmetric one, i.e. between all the nearest neighbors. For even N the next to the nearest neighbor pattern of stable exact synchronization is also possible. Patterns of exact synchronization in diffusively coupled type II systems depends on N being even or odd in a similar way as for the type I systems. For odd N the only pattern of exact synchronization is between all nearest neighbors and can be achieved either by increasing s or c. In the case of sigmoid coupling and small or zero time-lag all units can oscillate (for c > c0) only in exact synchrony independently of N. For larger s odd N implies only nearest neighbor exact synchronization, and for even N also the next to the nearest neighborÕ pattern of exact synchronization can be stable.
Acknowledgment This work is supported by the Serbian Ministry of Science contract No. 101443,101865,101225.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]
Rabinovitch A, Rogachevskii I. Threshold, excitability and isochrones in the Bonhoeffer–van der Pol system. Chaos 1999;9:880–6. Izhikevich EM. Neural excitability, spiking and bursting. Int J Bif Chaos 2000;10:1171–266. Terman D, Wang D. Global competition and local cooperation in a network of neural oscillators. Physica D 1995;81:148–76. FitzhHugh R. Mathematical models of threshold phenomena in the nerve membrane. Bull Math Biophys 1955;17:257–78. Nagumo J, Arimoto S, Yoshizava S. An active pulse transmission line simulating nerve axon. Proc IREE 1962;50:2061. Fox JJ, Jayaprakash C, Wang D, Campbell SR. Synchronization in relaxation oscillator networks with conduction delays. Neural Comput 2001;13:1003–21. Campbell SR, Wang D. Relaxation oscillators with time delay coupling. The Ohio State University CIS-Technical Report 1996, No. 47. Renvarez G. Synchronization in two neurons: results for two-component dynamical model with time-delayed inhibition. Physica D 1998;114:147–71. Ramana Reddy DV, Sen A, Johnston GL. Time delay induced death in coupled limit cycle oscillators. Phys Rev Lett 1998;80:5109–12. Ramana Reddy DV, Sen A, Johnston GL. Time delay effects on coupled limit cycle oscillators at Hopf bifurcation. Physica D 1999;129:15. Hoppensteadt FC, Izhikevich EM. Weakly connected neural networks. New York: Springer-Verlag; 1997. Izhikevich EM. Phase models with explicit time delay. Phys Rev E 1998;58:905–8. Schuster HG, Wagner P. Mutual entrainment of two limit cycle oscillators with time delayed coupling. Prog Theor Phys 1989;81:939–45. Nakamura Y, Tominaga F, Munakata T. Clustering behaviour of time-delayed nearest-neighbor coupled oscillators. Phys Rev E 1994;49:4849. Kim S, Park SH, Ryu CS. Multistability in coupled oscillator systems with time delay. Phys Rev Lett 1997;79:2911. Yeung MKS, Strogatz SH. Time delay in the Kuramoto model of coupled oscillators. Phys Rev Lett 1999;82:648–51. Grossberg S. Nonlinear difference-differential equations in prediction and learning theory. Proc Natl Acad Sci USA 1967;58:1329–34. Gopalsamy K, Leung I. Delay induced periodicity in a neural network of excitation and inhibition. Physica D 1996;89:395. Ye H, Michel A, Wang K. Qualitative analysis of Cohen–Grossberg neural network with multiple delays. Phys Rev E 1995;51:2611. Olier J, Belair J. Bifurcation, stability and monotonicity properties of a delayed neural network model. Physica D 1997;102:349–63. Shayer LP, Campbell SA. Stability, bifurcation, and multistability in a system of two coupled neurons with multiple time delays. SIAM J Appl Math 2000;61:673–700. Ncube I, Campbell SA, Wu J. Change in criticality of synchronous Hopf bifurcation in a multiply delayed neural system. Fields Ins Commun 2002;36:1–15. Buric´ N, Todorovic´ D. Dynamics of FitzHugh–Nagumo excitable systems with delayed coupling. Phys Rev E 2003;67:066222. Buric´ N, Grozdanovic´ I, Vasovic´ N. Excitable and oscillatory dynamics in an in-homogeneous chain of excitable systems with delayed coupling. Chaos, Solitons & Fractals 2004;22:731–40.
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
1233
[25] Murray JD. Mathematical biology. New York: Springer; 1990. [26] Hu B, Zhou C. Phase synchronization in coupled nonidentical excitable systems and array-enhanced coherence resonance. Phys Rev E 2000;61:R1001. [27] Pikovsky A, Rosenblum M, Kurths J. Synchronization. A universal concept in nonlinear sciences. Cambridge: Cambridge University Press; 2001. [28] Hale J, Lunel SV. Introduction to functional differential equations. New York: Springer-Verlag; 1993. [29] Pyragas K. Synchronization of coupled time-delay systems: analytical estimations. Phys Rev E 1998;58:3067–71.
Type I vs. type II excitable systems with delayed coupling Nikola Buric´
a,*
, Ines Grozdanovic´ b, Nebojsˇa Vasovic´
b
a
b
Department of Physics and Mathematics, Faculty of Pharmacy, University of Beograd, Vojvode Stepe 450, Beograd, Serbia and Montenegro Department of Applied Mathematics, Faculty of Mining and Geology, P.O. Box 162, Beograd, Serbia and Montenegro Accepted 8 June 2004
Abstract Excitable and oscillatory dynamics of delayed locally coupled type I and type II excitable systems is analyzed. Diffusive and sigmoid coupling have been considered. It is shown that the stability and the patterns of exactly synchronous oscillations depend on the type of excitability and the type of coupling. However, within the same class, characterized by the excitability type and the coupling the dynamics qualitatively depends only on whether the number of units is even or odd. 2004 Elsevier Ltd. All rights reserved.
1. Introduction Excitability is a common property of many physical and biological systems. Although there is no precise definition [1] the intuitive meaning is clear: A small perturbation from the single stable stationary state can result in a large and long lasting excursion away from the stationary state before the system is returned back asymptotically to equilibrium. Further more, as an external parameter is changed, the global attractor in the form of the stationary point bifurcates into a stable periodic orbit, and the excitability is replaced by the oscillatory dynamics. The are basically two qualitatively different types of excitability, fancifully called type I and type II, which are distinguished by different bifurcations leading to the oscillatory dynamics [2]. This is manifested in different properties of the frequency of the emergent oscillations. Our goal in this paper is to study the differences in the dynamics of relatively small chains of coupled type I or type II excitable systems. As the example of the type I excitability we shall use the Terman–Wang model [3] for the corresponding values of the parameters, and as the example of the type II excitable system we shall use the FitzHugh–Nagumo model [4,5]. Both systems originate as models or neuronal excitability, i.e. of the dynamics of the neuron transmembrane potential, but this interpretation of the models will not be of primary importance for us. We should stress that we shall analyze the systems such that each of the units in a chain is a system operating in the excitable, and not in the oscillatory, regime. The oscillations might appear only due to the coupling between the excitable units. Transfer of excitations between the neighboring units is not instantaneous, so that realistic models of the coupling should depend on the states of the neighboring units at different times. Thus we shall suppose that the function which
*
Corresponding author. E-mail address: [email protected] (N. Buric´).
0960-0779/$ - see front matter 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2004.06.033
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
1222
describes the influence of the ith unit on the jth unit at the time t depends on the state of the ith unit at some earlier time t s. In fact we shall consider two types of coupling: the so-called diffusive, described by f(xj, xi) = xj(t) xi(t s) and the sigmoid type given by f ðxj ; xi Þ ¼ tan1 xi ðt sÞ. The effects of time-lags on the dynamics of coupled excitable systems are much less studied than for the systems of delayed coupled regular or chaotic oscillators. Let us illustrate typically studied systems by quoting only some representative references. Delayed coupled relaxation oscillators have been studied for example in [6–8]. In these studies, each unit is a relaxation oscillator, and the primary objective of the analysis are the periodic orbits that appear in the delayed coupled system. Singular approximation, or an approximate or numerical construction of the Poincare map, are used to analyze various types of synchronous or asynchronous oscillations. In [9,10], a network of N = 2 and N > 2 oscillators described by the equations of the normal form of the Hopf bifurcations with delayed diffusive coupling is studied. It is shown that the time-delay can lead to the stabilization of the trivial stationary point, which is interpreted as the amplitude death. Less directly related to our work is the analyzes of the influence of the time delay in the systems of coupled phase oscillators. The influence of time-delay in coupled phase oscillators was studied for example in [11,12], where it is shown that the time-delay cannot introduce significant changes into the dynamics of a class of such systems, unless the time-lag is of the order of several oscillation periods [12]. Independently of neuronal models, collective behavior of the phase coupled (phase) oscillators with time-delayed coupling, have been studied using the dynamical [13–15], or statistical [16] methods. Finally, the influence of time delay has been studied in the Cohen– Grasberger–Hopfield type of neural networks, as early as in 1967 [17]. More recent references are for example [18,19], and for small networks [20,21] (see also [22] and the references there in). Delayed coupled type II excitable, as opposed to oscillatory, systems have been recently studied in reference [23] (where a collection of references to the typical results for coupled oscillators is given). This was extended in [24] to chains of FitzHugh–Nagumo systems with small in-homogeneities. In these papers only the type II system, given by the FitzHugh–Nagumo equations, and the sigmoid delayed coupling is considered. In this paper we shall compare the characteristic features of delayed coupled type I and type II systems for different typical forms of the coupling, and consider in more details the properties of the patterns of synchronous oscillations that occur in these systems. The paper is organized as follows. In Section 2 the two basic models of type I and type II systems and the two types of coupling will be introduced. In Section 3 the local bifurcations of the relevant stationary solutions will be analyzed. Stable oscillatory dynamics will be studied in Section 4. Section 5 contains the summary of our main conclusions.
2. Two basic models In this section we introduce the class of models that will be analyzed. We shall consider one-dimensional chains of N excitable units with periodic boundary conditions, where only the nearest neighbors are coupled. Such a system is given by the equations of the following general form: x_ i ¼ X ðxi ; y i Þ þ Cðxi ; xiþ1 ðt sÞ; xi1 ðt sÞÞ; y_ i ¼ Y ðxi ; y i Þ; i ¼ 1; 2; . . . ; N ;
ð1Þ
with the following boundary conditions on Nth and 0th unit: ðx0 ; y 0 Þ ¼ ðxN ; y N Þ;
ðxN þ1 ; y N þ1 Þ ¼ ðx1 ; y 1 Þ:
ð2Þ
The field X(x, y), Y(x, y) describes the single isolated unit, and the function Cðxi ðtÞ; xii ðt sÞ; xiþ1 ðt sÞÞ Cðxi ðtÞ; xsi1 ðtÞ; xsiþ1 ðtÞÞ describes the time-delayed coupling between the nearest neighbors. As the examples of type I and type II excitability we shall take the Terman–Wang and FitzHugh–Nagumo models, respectively. In the first case, the single unit is described by X ðx; yÞ ¼ X TW ¼ x3 þ 3x y þ 2 þ I; Y ðx; yÞ ¼ Y TW ¼ ½að1 þ tanhðx=bÞÞ y;
ð3Þ
where we fix the parameters to the following typical values a = 3, b = 0.2, I = 0.4, = 1, so that the system has a stable stationary solution (x0, y0) (1.34577, 0), which would go through the saddle-node on the limit cycle bifurcation upon increasing the absolute value of the external parameter I. The resulting stable limit cycle has finite radius and the frequency is vanishingly small. This type of bifurcation characterizes type I excitability. For type II excitability we shall take the FitzHugh–Nagumo model in the form [25]
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
1223
Fig. 1. The nulclines and a typical excitable orbit for the type I system (). The parameters are a = 3, b = 0.2, I = 0.4 and the stationary state is at (x0, y0) (1.34577, 0).
X ðx; yÞ ¼ X FN ¼ x3 þ ða þ 1Þx2 ax y þ I; Y ðx; yÞ ¼ Y FN ¼ bx cy:
ð4Þ
where a, b and c are positive internal parameters, and I, which will be taken I = 0, is the external parameter. The stationary solution (0, 0) is stable for 4b < c(a 1)2, I = 0. The stationary solution loses stability through the Hopf bifurcation upon increasing I and a stable limit cycle with vanishingly small radius but finite frequency is born. This type of bifurcation characterizes type II excitability. Notice that x_ nullcline in both models have the same quibic shape, and the main difference between the models is in y_ nullcline, which is linear for the FitzHugh–Nagumo and sigmoid for the Terman–Wang. The nullclines and typical excitable orbits are presented in Figs. 1 and 2, for the Terman–Wang and the FitzHugh–Nagumo systems, respectively. Two types of coupling will be analyzed. The diffusive coupling, Cðxi ðtÞ; xsi1 ðtÞ; xsiþ1 ðtÞÞ ¼ cð2xi þ xsi1 þ xsiþ1 Þ
ð5Þ
and the sigmoid coupling in the following form: Cðxi ðtÞ; xsi1 ðtÞ; xsiþ1 ðtÞÞ ¼ ctan1 ðxsi1 x0 Þ þ ctan1 ðxsiþ1 x0 Þ;
ð6Þ
where (x0, y0) is the stable stationary solution. Thus, the coupling between the units in equilibrium is zero. As expected, the model of the coupling plays important role in the dynamics of the whole chain, but, nevertheless, some properties of the synchronization patterns are independent of the type of the coupling.
3. Bifurcations of stationary states In this section we present results of the analyses of local bifurcations due to time-lag of the stationary solution for a system of only two delayed-coupled excitable units. Only the stationary solution which corresponds to the excitable dynamics will be considered. Coordinates of the stationary solution do not depend on the coupling constant c or the time-lag s, but the stability could be effected by varying either c or s. We shall be interested in the dependence
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
1224
Fig. 2. The nulclines and a typical excitable orbit for the type II system (). The parameters are a = 0.25, b = c = 0.02, and the stationary state is at (x0, y0) = (0, 0).
of the bifurcation values of the time-lag s on the coupling constant c, for all intrinsic parameters of each unit fixed to some typical values, such that the considered stationary solution is stable when c = 0 and s = 0. The bifurcations occur for the values of the parameters and s such that there is a non-hyperbolic root, k, Re k = 0 of the characteristic equation. As we shall see the relevant characteristic equation has the same general form for all considered models and couplings, the differences between the models occurring in the explicit relations between the parameters and the coordinates of the considered stationary solutions. Linearization of the any of the considered models, around an arbitrary stationary solution x10 = x20 ” x0 and y10 = y20 ” y0 and substitution: xi(t) = Ai exp kt; yi = Bi exp kt results in the following characteristic equation: DðkÞ ½ðA þ kÞð þ kÞ þ B2 ½cð þ kÞ expðksÞ2 ¼ 0;
ð7Þ
which can be factored as D(k) = D1(k)D2(k) where D1;2 ðkÞ ¼ ½k A þ B c expðksÞð þ kÞ þ B ¼ 0:
ð8Þ
In the Terman–Wang case we have A ¼ 3 3x20 ;
B¼
a bcosh2 ðx0 =bÞ
ð9Þ
for the sigmoid coupling, and for the diffusive coupling A is replaced by A c. In the FitzHugh–Nagumo case we have A a;
B b;
¼c
ð10Þ
for the sigmoid coupling, and for the diffusive coupling A is replaced by A c, analogously as for the type I case. The difference between the systems and the couplings is further manifested in the characteristic equation when the explicit expressions for the respective stationary solutions (x0, y0) are substituted. Let us consider first the zero solutions of (7). This type of bifurcations occurs also when s = 0. Substitution k = 0 implies
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
A B c ¼ 0:
1225
ð11Þ
Next we consider pure imaginary solutions of D(k) = 0. If some additional non-degeneracy conditions are satisfied this corresponds to the Hopf bifurcation. Substitution k = ix where x is assumed real and positive, into D1 = 0 or D2 = 0, and separation of the imaginary and real parts of the resulting equations gives a biquadratic equation for x, x4 þ Mx2 þ N ¼ 0;
ð12Þ
where M ¼ A2 þ 2 2B c2 ;
N ¼ A2 2 þ B2 2AB c2 2
with the formal solutions: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2 ¼ ðM M 2 4NÞ=2;
ð13Þ
ð14Þ
which of the solutions x± are real depends on the model and its stationary solution. For the stationary solutions which represent the quiescent states of the units, we shall see that, only the pair of FitzHugh–Nagumo systems with the sigmoid or the diffusive coupling do have both type of solutions, i.e. both x2þ and x2 are real and positive. In the case of the type I systems with either types of coupling, the only real solution, for the corresponding unexcited stationary state, is x+. The critical values of the time-lag sc as a function of the parameters is obtained by inverting: tanðsc Þ ¼
xð2 þ x2 BÞ ; A2 Ax2 B
ð15Þ
where the upper or lower signs correspond to D1/2 = 0, and substituting the real and positive solutions for x of (12) in (15). The nature of the Hopf bifurcation for s = sc is determined by the dependence of the real part of k on s as it goes through sc. This is given by the sign of d Re k=ds for s = sc. Using the factorized characteristic equation, one obtains d Re k A2 þ 2 2B þ 2x2 c2 ¼ : ð16Þ ds s¼sc c2 ðx2 þ 2 Þ The previous expression is positive for x = x+, and the Hopf bifurcation increases the number of unstable direction as s is increased through sc. In the case there is real and positive x the expression (16) is negative for x = x, and the Hopf bifurcation decreases the number of unstable directions as s is increased through sc. As we have pointed out, the latter case occurs only for the FitzHugh–Nagumo systems. In what follows we denote by scj; , j = 0, 1, 2, . . . , different branches of the bifurcation curve given by (14) and (15). Let us now apply the formulas (11), (14) and (15) to obtained the bifurcations curves in the (c, s) plane for particular models and the relevant stationary solutions. Thus, we fix the internal parameters of each unit to such values that the stationary solution which correspond to the unexcited state is stable. The values of the coordinates of the stationary solution (x0,y0) do not depend on the coupling constant c. Internal parameters and the corresponding (x0, y0) are substituted into the formulas for x and sc. The result is the set of bifurcation curves scj; for the considered stationary solution. 3.1. Type I systems In the case of the Terman–Wang model, we take for both units = 1, a = 3, I = 0.4, b = 0.2. The point (x1, 0, y1, 0, x2, 0, y2, 0) = (1.34577, 0, 1.34577, 0) approximates the stationary solution up to the fifth decimal place. The bifurcation curves for the sigmoid and the diffusive coupling are presented in Fig. 3a and b. We point out that the solution x+ is positive only for such values of the coupling constant c for which the stationary solution is unstable when s = 0. There is a set of destabilizing Hopf bifurcation curves scj;þ in the (c, s) plane. The solution x is always negative, so there are no stabilizing Hopf bifurcations. Thus, the local effect of the time-delay on the stationary solution is only an increase of the number of unstable directions of the considered stationary solution. On the other hand, stable periodic orbits could appear due to bifurcations of other stationary solutions, which are stable for s = 0 and the considered c, or due to some global bifurcation. Thus, the existence of stable oscillatory solutions is not indicated by the local analysis of the considered stationary solution, and will be analyzed by numerical calculations, as is presented in the next section.
1226
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
Fig. 3. Hopf bifurcation curves scj;þ ðcÞ for the type I system and (a) sigmoid (b) diffusive coupling. The fixed values of the local parameters are as in Fig. 1.
3.2. Type II systems The Hopf bifurcation curves for two diffusively or sigmoid coupled FitzHugh–Nagumo, i.e. type II systems, are presented in Fig. 4. The values of the internal parameters are for both units a = 0.25, b = c = 0.02. The stationary solution is (x1, 0, y1, 0, x2, 0, y2, 0) = (0, 0, 0, 0) and, due to the polynomial form of the equations, it is found exactly. The stationary solution is stable for these values of the internal parameters and jcj 6 c0. The case of sigmoid coupled FitzHugh–Nagumo units was thoroughly discussed in [23,24]. Important difference, as compared with type I systems, is that even the instantaneous coupling of excitable type II units leads to a supercritical Hopf bifurcation at c0 = (a + c)/d, c < 1, which results in a stable periodic orbit. For the sigmoid coupling, the dynamics of the two units on the resulting periodic orbits is exactly synchronized, i.e. x1(t) = x2(t); y1(t) = y2(t). Further important characteristic feature of this case is the existence of real x solutions and the corresponding stabilizing Hopf bifurcations on scj; , j = 1, 2, . . . As was discussed in the cited references, this together with a global fold limit cycle bifurcation leads to the phenomenon of the oscillator death in a certain domain of the (c, s) plane. The diffusive coupling has similar effect on the stability of the stationary solution. Again, there is a supercritical Hopf bifurcation due to the instantaneous coupling (s = 0), at some c = c0 = (ac + c) < 0, and the stationary solution becomes unstable. However, the stable periodic orbit is born in the plane x1 = x2;y1 = y2 and corresponds to the out of phase oscillations. Non-zero time-lag can lead to either stabilizing or destabilizing Hopf bifurcations. In fact, both solu-
Fig. 4. Hopf bifurcation curves scj;þ ðcÞ for the type II system and (a) sigmoid (b) diffusive coupling. The fixed values of the local parameters as in Fig. 2.
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
1227
tions x+ and x can be real. Again, the oscillator death occurs in a certain domain of the (c, s) plane. Let us stress again that the phenomenon of oscillator death does not occur for coupled type I systems.
4. Oscillatory dynamics Our goal in this section is to describe the main types of oscillatory solutions and to study the stability of synchronization between different units. Finally, we shall compare the conclusions about type I systems with the corresponding behavior of type II systems. In order to describe the stable periodic solutions of N coupled units we have calculated numerically the mean frequencies of the projections of the particular stable periodic orbit on (xi, yi) plane. Stability of the exact synchronization is studied by calculation of the largest transversal Lyapunov exponent. Let us briefly recapitulate the definitions of the mean frequency and the transversal Lyapunov exponent for DDE. The mean frequency of the ith unit is introduced via the corresponding phase which is commonly defined by [26,27] ui ðtÞ ¼ 2p
t tki þ 2pk; tki
tkþ1 i
ð17Þ
where tki is the time of the kth crossing of xi(t) of some convenient threshold. If the intersection times are eventually equally distributed then the curve ui(t) is asymptotic to a straight line, i.e. limt!1 ui ðtÞ=xi t ¼ 1. The constant xi is the mean frequency, with the corresponding period Ti = 2p/xi, of the ith unit. The concept of the mean frequency enables us to talk about coherent dynamics even when the oscillations of each unit are not harmonic. The dynamics of two units i and j is coherent if the mean frequencies are equal, i.e. xi = xj Two coherent units i and j oscillate out-of-phase if tki ¼ tkj þ px1 . If tki ¼ tkj the two units oscillate in-phase. In order to study the stability of the exact synchronization between ith and jth unit we numerically calculate the dynamics: Di;j ðtÞ ¼ xi ðtÞ xj ðtÞ:
ð18Þ
A stable trivial stationary solution Di, j = 0 corresponds to the stable exactly synchronous solution xi(t) = xj(t). It is useful (see for example [28,29] to consider the dynamics of Di, j as a dynamical system on the phase space given by continR0 uous functions D defined on the interval [s, 0], with the norm kDk2 ¼ s D2 ðhÞ dh. Solutions for different initial functions in C[s, 0] generate a semi-flow on this phase space given by Dt(h) = D(t + h), t 2 R+, h 2 [s, 0]. The largest transversal Lyapunov exponent for the (i, j) exact synchronization is then defined in analogy with the finite-dimensional systems by 1 kDi;j;t k Li;j ¼ lim ln ; ð19Þ t!1 t kDi;j;0 k where the norm iDi,j;0i is small. (i, j) synchronization is locally attracting for some (c, s) if Li, j is negative for at least some sufficiently small Di,j,0. 4.1. Type I systems Chains of type I systems with delayed diffusive or sigmoid coupling have quite different global behavior. We shall first discuss the case of the diffusive coupling. The main results of our computations in this case are presented in Figs. 5 and 6. The parameters a, b and I are fixed to some typical values (for example, the same as in the previous section). Many initial functions are evolved, and the transients are discarded. There are no oscillatory solutions for any c and s = 0 or sufficiently small. However, stable oscillatory behavior appears for the time-lag greater then some scr. This scr and the oscillatory solutions are in no way related to the Hopf bifurcations dynamics of the considered stationary solution which was studied in the previous section. In fact, oscillatory solutions are possible only for c such that the stationary solution is unstable and scr is larger then s0,+. Typical dependence of the mean frequency x1, of the first component x1 on the stable limit cycles, on the time-lag, for N = 5 and N = 6, is illustrated in Figs. 5a and b. In fact, we found that the mean frequencies of all components are equal, so that the oscillations are coherent for any coupling c and time-lag s. Thus, Figs. 5 represent the frequency of any of the components. Furthermore, the typical dependence of the common frequency on the time-lag is approximately independent of N. The general conclusion is that for c < 0 the mean frequency decreases as s is increased. For c > 0 and s larger than some critical value scr the frequency first strongly depends on s in a small interval beyond scr.
1228
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
Fig. 5. Mean frequency xi as a function of s for diffusively coupled type I systems: (a) N = 5 and (b) N = 6. The fixed values of c are (1) c = 1 (boxes), (2) c = 4 (crosses), (3) c = 4 (triangles). The local parameters are as in Fig. 1.
Fig. 6. Coherent oscillations for diffusively coupled N = 6 type I systems: (a) out of phase (xi, xi + 1) and exactly synchronous (xi, xi + 2) for c = 1, s = 3; (b) convergence of the computations of the transversal Lyapunov exponents for c = 4, s = 20; (c) oscillations for typical parameters; (d) transient dynamics for large s = 50.
Further away from scr, the frequency oscillates with s in a quite narrow band. For small variable c, and fixed s, the frequency increases with c, and then weakly depends on the larger c. Let us stress that the same behavior is observed for any N as large as N = 50. The patterns of stable synchronization depend on N being odd or even. For any odd N the nearest neighbors always oscillate in exact synchrony, which is the simplest spatially uniform pattern. This is also the only patters that we have observed for any even N and c < 0. However, for even N, and c > 0, next to the nearest neighbors stable exactly synchronous oscillations are possible, i and i + 1 unit could oscillate out of phase, as it is illustrated in Fig. 6a. The out of
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
1229
phase dynamics between the nearest neighbors was observed only in a narrow interval of parameters c and s scr(c). In general exact synchronization of the nearest neighbors is much slower than the synchronization between the next to the nearest units. This is clearly shown by calculations of the transverse Lyapunov exponent Li;j ¼ limt!1 Li;j ðtÞ, illustrated in Fig. 6b. Let us point out that, for larger time-lags, say s = 50 like in Fig. 6d, the approach to the periodic oscillations on the attractor, could involve quite long transient dynamics, during which a sequence of almost synchronous bursts of oscillations alternate periods of stationary behavior. One could contemplate a source of a bursting behavior, which is executed by appropriate variations of the time-lag in the transmission of interaction through the junctions (synapses) between the units. The behavior of chains with type I units and delayed sigmoid coupling like in (6) is trivial. The symmetric stationary solution, that corresponds to the equilibrium, remains stable for any s, as long as the coupling c is such that the state is stable for s = 0. Larger values of c destabilize this state, new stable stationary states are introduced, but no oscillations occur. 4.2. Type II systems Instantaneous diffusive coupling can induce oscillations in a chain of type II excitable systems via Hopf bifurcation for the coupling parameter c = c0 < 0. There are no oscillations for c > 0. The resulting oscillations, for c < c0 < 0, of each unit are coherent. The transition from regular monochromatic oscillations of each unit, for c c0 to relaxation oscillations occurs in a very small interval of c < c0. Nevertheless the mean frequency of different units, as defined in (17), is the same for any c. Let us point out that if the parameter c is such that the units are excitable, i.e. the only attractor is the stable stationary solution, and c is far away from c0, then the time-lag has no effect on the qualitative properties of the dynamics. The stationary solution remains stable for any s. There are two distinct domains of the coupling c for which the dynamics is oscillatory. Firstly, there is a small interval near c0 for which the instantaneously coupled units are excitable but the time-lag does introduce stable oscillations, as it is indicated by the bifurcations curves in Fig. 4. The other situations occurs for jcj > jc0j, i.e. when the instantaneously coupled units oscillate, and when further Hopf bifurcations are possible due to non-zero time-lag. The results of our numerical computations are presented in Figs. 7–9. Fig. 7a and b, illustrate the dependence of the mean frequency on the time-lag for the values of c very close to c0, and for some typical value further away from c0. The Hopf bifurcation of the stationary solution at c0 (for s = 0) is responsible for the appearance of the oscillatory dynamics, so that the frequency near the bifurcation has a finite non-zero value. The dependence of xi(s) for fixed c is piecevise monotonically decreasing. The bifurcation values of s, corresponding to the jumps in xi(s), also correspond to changes of the patterns of synchronization between the units, i.e. to the threshold of the nearest neighbor exact synchronization. Possible patterns are illustrated in Figs. 8 and 9. The jumps do not occur if the fixed value of c is very close to c0. Again, the chains with odd or even N behave differently. For odd N and jcj > jc0jclose to, or moderately away from c0, and for small s, there is no stable exact synchronization between any two units. Larger s leads to stable exact synchronization. However, the only pattern of the stable exact synchronization is the most symmetric one, i.e. between all of the units. The synchronization threshold of the time-lag s decreases with increasing the absolute value of c < c0. For example, for c = 0.1 the threshold is s = 22 and for c = 0.3 the threshold is s = 9. Sufficiently large c induces exactly
Fig. 7. Mean frequency xi as a function of s for diffusively coupled type II systems: (a) N = 5 and (b) N = 6. The fixed values of c are (1) c = 0.1 (boxes), (2) c = 0.3 (crosses), (3) c = 0.056 (triangles). The local parameters are as in Fig. 2.
1230
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
Fig. 8. Illustrates transition from asynchronous coherent to exactly synchronous (i, i + 1) pattern as s is increased, for diffusively coupled N = 5 type II systems (). Frame 8f shows the transversal Lyapunov exponent for s = 20 (at the synchronization threshold) and s = 22 (slightly above the threshold) The local parameters are as in Fig. 1, c = 0.1.
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
1231
Fig. 9. Illustrates bifurcations of the synchronization patterns for diffusively coupled N = 6 type II systems (): (a) increasing c for fixed s = 20 leads to (xi, xi + 1) exact synchronization; (b) illustrates (xi, xi + 2) exact synchronization for s = 10 and c = 0.1.
synchronous oscillations with arbitrary small s. The synchronization threshold is clearly indicated by the transversal Lyapunov exponent. Chains with even N display two different patterns of stable exact synchronization. For the coupling c smaller than c0 but very close, the next to the nearest neighbors exactly synchronous oscillations are stable. The nearest neighbors oscillate out of phase. This pattern persists for any time-lag s. On the other hand, if the absolute value of c < c0 is larger, (i, i + 2) pattern is substituted by the simplest pattern with all units exactly synchronized. It is interesting to observe that for moderate values of s there could coexist locally stable attractors that correspond to different patterns of exact synchronization. For example, for N = 6, c = 0.1 and s = 15 we have found that the (i, i + 2) pattern is locally stable, but also the pattern with x1 = x4 and x6 = x3 = x2 has a non-zero domain of stability. Larger time-lag again leads to globally stable (i, i + 1) pattern of synchronization. Finally, we would like to make just a few comments on the oscillations and the patterns of synchronization in chains with sigmoid delayed coupling. Oscillations in chains of type II excitable units with such coupling have been analyzed before in [23,24]. Here we would like just to point out that the main differences compared with the diffusive coupling is that for zero time-lag sufficiently strong sigmoid coupling leads to stable exactly synchronous oscillations of all units. Thus we have to consider two different domains of coupling in analyzing the influence of the time-lag. If the instantaneously coupled units have excitable, and not oscillatory dynamics, large time delay can lead to a global fold limit cycle bifurcation, with a bi-stability as the result. The stationary state and a limit cycle are the two attractors. The oscillations between the neighbors are out of phase but the next to the nearest units oscillate in exact synchrony. Sufficiently strong instantaneous sigmoid coupling leads to exactly synchronous oscillations of all type II units, but already relatively small time-lags introduce the next to the nearest neighbor pattern for even N. As we have seen, the diffusive coupling leads to the opposite situation. Furthermore, the oscillations introduced by instantaneous coupling of type II systems, can be turned of by the time-delay. There is no such phenomenon in the chains of type I systems with either of the considered types of coupling.
5. Summary and discussion We have studied dynamics of chains of bidirectionally and locally coupled excitable systems. Sigmoid and diffusive coupling, with explicitly included time-delay, have been considered. Our main goal was to compare the behavior of type I vs. type II excitable systems with respect to the following two questions. The first one was the change of stability, with respect to coupling strength and time-delay, of the state that represent the attractor of the excitable dynamics. The second question was concerned with possible spatial patterns of exact synchronization between the oscillations of each unit that originate only through coupling. Let us stress that the isolated units are always only excitable and any oscillatory behavior is introduced by coupling. Concerning the first question, we have concluded that even the instantaneous coupling of either type can destabilize the stationary solution of the excitable behavior. However, for type II systems the destabilization occurs via Hopf bifurcation and a small stable limit cycle appears in a neighborhood of the considered stationary solution. On the other hand, the destabilization of the type I stationary state by the instantaneous interaction just introduces new stable stationary solution. In correspondence with this difference between the two types of excitable systems, the time-delay has
1232
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
different effect on the dynamics near the stationary solution. For type I systems the time-delay can only increase the number of unstable directions via supercritical Hopf bifurcation. For type II systems the time-delay can also decrease the number of unstable directions via inverse subcritical Hopf bifurcation. In particular, the time-delay can induce the oscillators death only in the chains of type II systems. Concerning the second question, we first point out that oscillatory solutions in chains of type I excitable systems are observed only for non-zero time-lag and diffusive coupling. In this case and for odd N, the only type of spatial pattern of exact synchronization is the most symmetric one, i.e. between all the nearest neighbors. For even N the next to the nearest neighbor pattern of stable exact synchronization is also possible. Patterns of exact synchronization in diffusively coupled type II systems depends on N being even or odd in a similar way as for the type I systems. For odd N the only pattern of exact synchronization is between all nearest neighbors and can be achieved either by increasing s or c. In the case of sigmoid coupling and small or zero time-lag all units can oscillate (for c > c0) only in exact synchrony independently of N. For larger s odd N implies only nearest neighbor exact synchronization, and for even N also the next to the nearest neighborÕ pattern of exact synchronization can be stable.
Acknowledgment This work is supported by the Serbian Ministry of Science contract No. 101443,101865,101225.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]
Rabinovitch A, Rogachevskii I. Threshold, excitability and isochrones in the Bonhoeffer–van der Pol system. Chaos 1999;9:880–6. Izhikevich EM. Neural excitability, spiking and bursting. Int J Bif Chaos 2000;10:1171–266. Terman D, Wang D. Global competition and local cooperation in a network of neural oscillators. Physica D 1995;81:148–76. FitzhHugh R. Mathematical models of threshold phenomena in the nerve membrane. Bull Math Biophys 1955;17:257–78. Nagumo J, Arimoto S, Yoshizava S. An active pulse transmission line simulating nerve axon. Proc IREE 1962;50:2061. Fox JJ, Jayaprakash C, Wang D, Campbell SR. Synchronization in relaxation oscillator networks with conduction delays. Neural Comput 2001;13:1003–21. Campbell SR, Wang D. Relaxation oscillators with time delay coupling. The Ohio State University CIS-Technical Report 1996, No. 47. Renvarez G. Synchronization in two neurons: results for two-component dynamical model with time-delayed inhibition. Physica D 1998;114:147–71. Ramana Reddy DV, Sen A, Johnston GL. Time delay induced death in coupled limit cycle oscillators. Phys Rev Lett 1998;80:5109–12. Ramana Reddy DV, Sen A, Johnston GL. Time delay effects on coupled limit cycle oscillators at Hopf bifurcation. Physica D 1999;129:15. Hoppensteadt FC, Izhikevich EM. Weakly connected neural networks. New York: Springer-Verlag; 1997. Izhikevich EM. Phase models with explicit time delay. Phys Rev E 1998;58:905–8. Schuster HG, Wagner P. Mutual entrainment of two limit cycle oscillators with time delayed coupling. Prog Theor Phys 1989;81:939–45. Nakamura Y, Tominaga F, Munakata T. Clustering behaviour of time-delayed nearest-neighbor coupled oscillators. Phys Rev E 1994;49:4849. Kim S, Park SH, Ryu CS. Multistability in coupled oscillator systems with time delay. Phys Rev Lett 1997;79:2911. Yeung MKS, Strogatz SH. Time delay in the Kuramoto model of coupled oscillators. Phys Rev Lett 1999;82:648–51. Grossberg S. Nonlinear difference-differential equations in prediction and learning theory. Proc Natl Acad Sci USA 1967;58:1329–34. Gopalsamy K, Leung I. Delay induced periodicity in a neural network of excitation and inhibition. Physica D 1996;89:395. Ye H, Michel A, Wang K. Qualitative analysis of Cohen–Grossberg neural network with multiple delays. Phys Rev E 1995;51:2611. Olier J, Belair J. Bifurcation, stability and monotonicity properties of a delayed neural network model. Physica D 1997;102:349–63. Shayer LP, Campbell SA. Stability, bifurcation, and multistability in a system of two coupled neurons with multiple time delays. SIAM J Appl Math 2000;61:673–700. Ncube I, Campbell SA, Wu J. Change in criticality of synchronous Hopf bifurcation in a multiply delayed neural system. Fields Ins Commun 2002;36:1–15. Buric´ N, Todorovic´ D. Dynamics of FitzHugh–Nagumo excitable systems with delayed coupling. Phys Rev E 2003;67:066222. Buric´ N, Grozdanovic´ I, Vasovic´ N. Excitable and oscillatory dynamics in an in-homogeneous chain of excitable systems with delayed coupling. Chaos, Solitons & Fractals 2004;22:731–40.
N. Buric´ et al. / Chaos, Solitons and Fractals 23 (2005) 1221–1233
1233
[25] Murray JD. Mathematical biology. New York: Springer; 1990. [26] Hu B, Zhou C. Phase synchronization in coupled nonidentical excitable systems and array-enhanced coherence resonance. Phys Rev E 2000;61:R1001. [27] Pikovsky A, Rosenblum M, Kurths J. Synchronization. A universal concept in nonlinear sciences. Cambridge: Cambridge University Press; 2001. [28] Hale J, Lunel SV. Introduction to functional differential equations. New York: Springer-Verlag; 1993. [29] Pyragas K. Synchronization of coupled time-delay systems: analytical estimations. Phys Rev E 1998;58:3067–71.