Synchronization of networks of oscillators with distributed delay coupling Y. N. Kyrychko, K. B. Blyuss, and E. Schöll Citation: Chaos: An Interdisciplinary Journal of Nonlinear Science 24, 043117 (2014); doi: 10.1063/1.4898771 View online: http://dx.doi.org/10.1063/1.4898771 View Table of Contents: http://scitation.aip.org/content/aip/journal/chaos/24/4?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Synchronization and array-enhanced resonances in delayed coupled neuronal network with channel noise Chaos 24, 033131 (2014); 10.1063/1.4894463 Phase locking of two limit cycle oscillators with delay coupling Chaos 24, 023123 (2014); 10.1063/1.4881837 Design of coupling for synchronization in time-delayed systems Chaos 22, 033111 (2012); 10.1063/1.4731797 Delay time modulation induced oscillating synchronization and intermittent anticipatory/lag and complete synchronizations in time-delay nonlinear dynamical systems Chaos 17, 013112 (2007); 10.1063/1.2437651 Synchronizing weighted complex networks Chaos 16, 015106 (2006); 10.1063/1.2180467
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 157.211.3.15 On: Mon, 12 Jan 2015 02:35:02
CHAOS 24, 043117 (2014)
Synchronization of networks of oscillators with distributed delay coupling €ll2 Y. N. Kyrychko,1,a) K. B. Blyuss,1 and E. Scho 1
Department of Mathematics, University of Sussex, Falmer, Brighton BN1 9QH, United Kingdom Institut f€ ur Theoretische Physik, Technische Universit€ at Berlin, Hardenbergstr. 36, 10623 Berlin, Germany
2
(Received 16 June 2014; accepted 9 October 2014; published online 24 October 2014) This paper studies the stability of synchronized states in networks, where couplings between nodes are characterized by some distributed time delay, and develops a generalized master stability function approach. Using a generic example of Stuart-Landau oscillators, it is shown how the stability of synchronized solutions in networks with distributed delay coupling can be determined through a semi-analytic computation of Floquet exponents. The analysis of stability of fully synchronized and of cluster or splay states is illustrated for several practically important choices of C 2014 AIP Publishing LLC. delay distributions and network topologies. V [http://dx.doi.org/10.1063/1.4898771] Many coupled real-life systems exhibit complex dynamical regimes, including the emergence of coherent dynamics known as synchronization. An important practical issue is understanding for what types of networks the synchronization can occur and is stable, which determines whether or not it can be observed experimentally. Quite often, connections between different elements of the network are non-instantaneous, and one has to explicitly include time delays in these connections in the analysis of synchronization. In this paper, we show how one can analyze stability of different types of synchronization in networks with distributed delay coupling, which represents a realistic situation in which the time delays themselves may not be constant or fixed. Our results suggest that the parameter space and network topologies for which a stable synchronization occurs are affected not only by the mean time delay but also by the width of the delay distribution.
I. INTRODUCTION
In the last decade, there has been a substantial growth of research interest in the dynamics of coupled systems, ranging from a few elements to large networks.1–4 From a perspective of potential applications, one of the central research questions for such systems is the emergence and stability of different types of collective dynamics and, in particular, synchronization.5,6 In order to study the stability of synchronization, Pecora and Carroll7 put forward a master stability function (MSF) approach that allows one to separate the local dynamics of individual nodes from the network topology, which is achieved by a proper diagonalization of the matrix representing the full network dynamics. An important aspect of network dynamics concerns the fact that interactions between nodes are often noninstantaneous due to a finite speed of signal propagation. In order to account for this feature, one has to explicitly include time delays in the analysis, and this is known to have a a)
Electronic mail: [email protected]
1054-1500/2014/24(4)/043117/10/$30.00
significant impact on network behaviour.8–20 Some work has been done on the extension of the master stability function to networks with time-delayed coupling10,21–26 but so far it has only included single constant time delays. Hunt et al.27,28 have considered complex networks of linearly coupled nodes with time delays and noise and established conditions for synchronizability of such networks. At the same time, in many realistic systems, time delays may be not constant and may either vary depending on the values of system variables, or just not be explicitly known,29–31 and in these cases, the standard methodology of considering one or several constant time delays is not sufficient. To describe such situations mathematically, one can use the formalism of distributed time delays, where the time delay is represented through an integral kernel describing a particular delay distribution.32–35 Distributed time delays have been successfully used to describe situations when only an approximate value of time delay is known in engineering experiments,35–37 for modelling distributions of waiting times in epidemiological models,38 maturation periods in population and ecological models,39,40 as well as in models of traffic dynamics,41 and neural systems.42–44 In a recent paper, Morarescu et al.45 have looked at systems with gamma-distributed delay coupling and analyzed the stability of the synchronized equilibria (steady states), which are not affected by the coupling. In this paper, we consider networks of coupled identical dynamical systems with linear distributed delay coupling as represented by some integral kernel. Linearization near the synchronization manifold yields a variational equation, which after block-diagonalization of the coupling matrix reduces to a single complex-valued integro-differential equation, whose maximum Floquet or Lyapunov exponent gives the master stability function in terms of the system parameters. Using the example of a network of Stuart-Landau oscillators, we make further analytical progress by showing how the computation of the master stability function can be reduced to finding solutions of a single transcendental characteristic equation. For particular choices of the delay distribution kernel, all terms in this equation can be found in a closed analytical form in terms of coupling and system parameters. The methodology we develop is quite generic and
24, 043117-1
C 2014 AIP Publishing LLC V
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 157.211.3.15 On: Mon, 12 Jan 2015 02:35:02
043117-2
€ll Kyrychko, Blyuss, and Scho
Chaos 24, 043117 (2014)
can be applied to analyze the stability of fully synchronized and of cluster states in any systems with distributed delay coupling. The outline of this paper is as follows. In Sec. II, we develop a general master stability function formalism for networks of identical coupled elements with distributed delay coupling. Section III illustrates how this general framework can be made more explicit by means of an amplitude-phase representation for the example of coupled Stuart-Landau oscillators. Section IV contains the results of analytical and numerical computations of the master stability function for coupled Stuart-Landau oscillators with uniform and gamma distributed delays and different coupling topologies. In Sec. V, the methodology of the master stability function for systems with distributed delay coupling is extended to the analysis of cluster states representing generalized synchronization. The paper concludes with discussion of results in Sec. VI. II. MASTER STABILITY FUNCTION
We consider a network of N identical dynamical systems with distributed delay coupling x_ k ¼ Fk ½xk ðtÞ þ r
N X
ð1 Gkj Hkj
0
0
0
gðt Þxj ðt t Þdt xk ðtÞ ;
are identical for different nodes, i.e., Hkj ¼ H for k; j ¼ 1; …N. The last two assumptions ensure that all nodes received the same input. With those assumptions, the model (1) simplifies to ð1 N X 0 0 0 Gkj H gðt Þxj ðt t Þdt xk ðtÞ ; x_ k ¼ F½xk ðtÞ þ r k ¼ 1; …; N:
k ¼ 1; …; N;
(1)
(3)
Let us now consider the synchronization manifold, defined as xk ðtÞ xs ðtÞ for all k ¼ 1; …N. The dynamics on the synchronization manifold is given by ð1 0 0 0 gðt Þxs ðt t Þdt xs ðtÞ ; (4) x_ s ¼ F½xs ðtÞ þ rlH 0
where l is the row sum, Eq. (2). It is worth noting that the dynamics of synchronization manifold is independent of the coupling strength in the case of zero row sum l ¼ 0 or instantaneous coupling gðuÞ ¼ dðuÞ, and in all other cases it will also depend on the coupling strength r. To study the stability of the synchronization manifold, we linearize the full system (3) near xk ðtÞ ¼ xs ðtÞ, which yields with xk ðtÞ ¼ xs ðtÞ þ nk ðtÞ,
0
j¼1
0
j¼1
n_ k ðtÞ ¼ J0 ðtÞnk ðtÞ þ r
N X j¼1
m
where xk 2 R is the state vector of each system, Fk ½ðxk ðtÞ describes the local dynamics of the node k, r 2 C is the coupling strength, G is an N N matrix describing the coupling topology and the strength of each link in the network, Hkj is an m m matrix representing the coupling scheme, and gðÞ is a delay distribution kernel, satisfying ð1 gðuÞdu ¼ 1: gðuÞ 0; 0
For gðuÞ ¼ dðuÞ, one recovers instantaneous coupling ðxj xk Þ; for gðuÞ ¼ dðu sÞ, the coupling takes the form of a discrete time delay: ½xj ðt sÞ xk ðtÞ. The matrix G, whose non-zero entries Gij define a relative strength of a link between nodes j and i, is related to the adjacency matrix of the network topology, whose elements are zero or one depending on the existence of the respective link. Next, we discuss a number of assumptions regarding the individual dynamical systems and their interactions in the network. All nodes are assumed to be identical, so that their dynamics is described by a single function Fk ½ ¼ F½. The time delays in the propagation of signals from different nodes are also taken to be the same and described by the delay distribution kernel g(u). Since we are interested in synchronization and its stability, we will further assume that the matrix G has a constant row sum N X
Gkj ¼ l:
(2)
ð1 Gkj H 0
k ¼ 1; …; N;
(5)
where J0 ðtÞ ¼ DF½xs ðtÞ is the Jacobian of F. Introducing a vector n ¼ ðn1 ; n2 ; …; nN ÞT , the above equation can be rewritten as ð1 _n ¼ IN ½J0 ðtÞ rlH n þ r½G H gðt0 Þnðt t0 Þdt0 ; 0
(6) where denotes the Kronecker product, and IN is the N N unity matrix. We will assume that matrix G is diagonalizable, i.e., there exists a unitary transformation U such that UGU1 ¼ diagðl; 1 ; 2 ; …; N1 Þ: Here, the first eigenvalue of the matrix G, which is equal to the row sum 0 l, corresponds to the eigenvector ð1; 1; …; 1Þ describing the dynamics on the synchronization manifold. Hence, this eigenvalue is called the longitudinal eigenvalue of G, whereas all other eigenvalues are called transverse eigenvalues as they correspond to directions transverse to the synchronization manifold. To make further analytical progress, it is instructive to diagonalize the coupling matrix G, which results in a blockdiagonalized variational equation
j¼1
This condition is necessary for the existence of the synchronized solution. Finally, we assume that the coupling schemes
gðt0 Þnj ðt t0 Þdt0 nk ðtÞ ;
f_ k ðtÞ ¼ ½J0 ðtÞ rlHfk ðtÞ þ rk H
ð1
gðt0 Þfk ðt t0 Þdt0 ;
0
k ¼ 0; …; N 1;
(7)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 157.211.3.15 On: Mon, 12 Jan 2015 02:35:02
043117-3
€ll Kyrychko, Blyuss, and Scho
Chaos 24, 043117 (2014)
where k is the k-th eigenvalue of the coupling matrix G. It is worth noting that in this system of equations, the Jacobian J0 ðtÞ and H are the same for all k, and one can consider the variational equation in dependence on the complex parameter w þ ib ð1 _fðtÞ ¼ ½J0 ðtÞ rlHfðtÞ þ ðw þ ibÞH gðt0 Þfðt t0 Þdt0 : 0
(8) Computation of the maximum Lyapunov exponent (or real part of the Floquet exponent) for this equation yields the master stability function
Introducing amplitude and phase variables as rk ¼ jzk j and /k ¼ argðzk Þ, Eq. (9) can be recast as a system of 2N real equations ð1 N X gðt0 Þrj ðt t0 Þcos r_ k ¼ k rk2 ðtÞ rk ðtÞ þ K Gkj 0
j¼1
i /j ðt t0 Þ /k ðtÞ þ h dt0 rk ðtÞcos h ; h
ð1 N X rj ðt t 0 Þ sin Gkj gðt0 Þ rk ðtÞ 0 j¼1 0 0 ½/j ðt t Þ /k ðtÞ þ hdt sin h ; k ¼ 1; …; N:
/_ k ¼ x crk2 þ K
Re½Kmax ðw; b; rÞ
(10)
and the analysis of stability of the synchronization manifold reduces to checking that for all transverse eigenvalues of the coupling matrix rk ¼ w þ ib; k ¼ 1; …; N 1, one has Re[Kmax ðw; b; rÞ < 0. A positive value of the Lyapunov exponent (real part of the Floquet exponent) associated with the longitudinal eigenvalue 0 indicates that the synchronized solution is chaotic (unstable). The significant advantage of this approach lies in the separation of the actual dynamics from the topology of the coupling, as the computation of Kmax ðw; b; rÞ can be done once for Eq. (8) independently of the coupling matrix. It is worth noting that although the computation of the master stability function can be performed independently of the coupling topology, it has to be done separately for each value of the coupling strength r. The reason for this is the fact that the row sum l is, in general, non-zero, and also the synchronized state xs ðtÞ itself depends on the coupling strength due to the delayed nature of the coupling, as follows from Eq. (4). The above approach can be generalized in a straightforward manner to more complex forms of the coupling function H. III. COUPLED STUART-LANDAU OSCILLATORS
In order to illustrate the analysis of the stability of synchronization, we consider an array of N identical diffusively coupled Stuart-Landau oscillators (k ¼ 1; …; N) z_ k ðtÞ ¼ ðk þ ixÞzk ðtÞ ð1 þ icÞjzk ðtÞj2 zk ðtÞ ð1 N X 0 0 0 þr Gkj gðt Þzj ðt t Þdt zk ðtÞ ; j¼1
(9)
0
K 2 Rþ ;
rk ðtÞ ¼ r0 ;
/k ¼ Xt;
k ¼ 1; …; N;
(11)
with the common radius r0 and the common frequency X of oscillations being determined by the solutions of the following equations: r02 ¼ k þ Kl½Fc ðX; hÞ cos h;
(12)
X ¼ x cr02 þ Kl½Fs ðX; hÞ sin h; where we have introduced the auxiliary quantities ð1 gðt0 Þ cosðat0 þ bÞdt0 ; Fc ða; bÞ ¼ 0 ð1 gðt0 Þ sinðat0 þ bÞdt0 : Fs ða; bÞ ¼
(13)
0
To analyse the stability of the fully synchronized solution (11), we use a slightly modified ansatz for small perturbations around this solution, namely, rk ðtÞ ¼ r0 ½1 þ drk ðtÞ;
/k ðtÞ ¼ Xt þ d/k ðtÞ;
k ¼ 1; …; N; which yields a variational equation of the form n_ k ¼ ðJ0 KlR0 Þnk þ K
N X
ð1 Gkj
k ¼ 1; …; N;
(14)
where nk ¼ ðdrk ; d/k ÞT , and the matrices J0 ; R0 and R1 are given by
h 2 R;
where K and h are the strength and the phase of coupling, respectively. Such complex-valued couplings, which are equivalent to a rotational matrix H if zk is written in terms of real and imaginary part, have been shown to be important in overcoming the odd-number limitation of time-delayed feedback control,46 in controlling amplitude death,47,48 and in anticipating chaos synchronization.49
gðt0 ÞR1 ðt0 Þnj ðt t0 Þdt0 ;
0
j¼1
which is a prototype of dynamics near a supercritical Hopf bifurcation. Here, zk 2 C, k, x 6¼ 0, and c are real constants, where x is the intrinsic oscillation frequency, and r ¼ Keih ;
The fully synchronized solution of the system (10) is given by
J0 ¼ R0 ¼ R1 ðt0 Þ ¼
2r02
0
2cr02
0
! ;
Fc ðX; hÞ Fs ðX; hÞ
; Fs ðX; hÞ Fc ðX; hÞ cosðh Xt0 Þ sinðh Xt0 Þ sinðh Xt0 Þ
cosðh Xt0 Þ
:
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 157.211.3.15 On: Mon, 12 Jan 2015 02:35:02
043117-4
€ll Kyrychko, Blyuss, and Scho
Chaos 24, 043117 (2014)
Fc ða; bÞ ¼ FLc ða; b; 0Þ;
Equivalently, the above system can be rewritten as n_ ¼ IN ðJ0 KlR0 Þn ð1 gðt0 Þ½G R1 ðt0 Þ n ðt t0 Þdt0 ; þK
(15)
Fs ða; bÞ ¼ FLs ða; b; 0Þ:
Equation (18) generalises an earlier result of Choe et al.,21,22 which considered the master stability function for the case of a single discrete time delay.
0
with the 2N-dimensional vector n ¼ ðn1 ; …; nN ÞT . Diagonalizing the matrix G results in a block-diagonalized variational equation ð1 gðt0 ÞR1 ðt0 Þfk ðt t0 Þdt0 : f_ k ðtÞ ¼ ðJ0 KlR0 Þfk ðtÞ þ Kk 0
(16) Since in these equations, none of the matrices depends explicitly on time t, it is possible to find Floquet exponents as the eigenvalues K 2 C of the following characteristic equation ð1 gðt0 ÞR1 ðt0 ÞeKt0 dt0 ¼ 0; det J0 KlR0 KI2 þ Kk 0
(17) where I2 is a 2 2 identity matrix. More explicitly, this transcendental equation has the form K
2
þ 2½r02 þ KðlFc ðX; hÞ k FLc ðX; h; KÞÞK þ 2r02 K½lFc ðX; hÞ k FLc ðX; h; KÞ þ2cr02 K½lFs ðX; hÞ k FLs ðX; h; KÞ þ K 2 ½lFc ðX; hÞ k FLc ðX; h; KÞ2 þK 2 ½lFs ðX; hÞ k FLs ðX; h; KÞ2 ¼ 0;
(18)
where by analogy with (13), we have introduced the quantities ð1 0 gðt0 Þ cosðat0 þ bÞezt dt0 ; FLc ða; b; zÞ ¼ 0 ð1 0 gðt0 Þ sinðat0 þ bÞezt dt0 ; FLs ða; b; zÞ ¼
IV. COMPUTATION OF MASTER STABILITY FUNCTION
In order to illustrate how the master stability function can be used for the analysis of stability of a synchronization manifold, we first have to specify a particular topology of the network G, as well as the distributed delay kernel g(u). We will consider a number of different coupling topologies, each of which is characterized by a particular spectrum of eigenvalues fk g, as shown in Fig. 1. Assuming that the constant row sum l is different from zero, we normalize matrices G by dividing all elements by l, which results in a row sum equal to unity. The coupling matrices G and their eigenvalues are then given by 0 1 0 1 0 0 0 B0 0 1 0 0C B C B. . . . .. C uni 2pik=N B . . . . ; G ¼ B . . . . . C C; k ¼ e B C @0 0 0 0 1A 1
0 0
0 (19)
for uni-directional ring coupling, 0 1 0 1 0 0 1 B C B1 0 1 0 0C B C B0 1 0 1 0C B C 1 Gbi ¼ B .. C .. .. .. .. C; 2B B . . . . . C B C B0 0 0 0 1C @ A 1 0 0 0 0
k ¼ cos
2pk ; N
k ¼ 0; …; N 1;
(20)
for bi-directional ring coupling, and
0
where the superscript L refers to these integrals representing the Laplace transform of modified kernels gðsÞ cosðas þ bÞ and gðsÞ sinðas þ bÞ. Comparing these expressions to (13) yields the following relations:
0
k ¼ 0; …; N 1;
Gall ij
8
1 (known as strong delay kernel in the case p ¼ 2), the biggest influence on the coupling at any moment of time t is from the values of zk at t ðp 1Þ=a. The delay distribution (30) has the mean time delay ð1 p (31) ugðuÞdu ¼ ; sm ¼ a 0 and the variance
ð1mÞþN mod N cm x k ;
Var ¼
m¼1
k ¼ 0; …; N 1;
ð1
gðuÞ ¼
k ¼ c1 þ cN xk þ cN1 x2k þ þ c2 xN1 k ¼
ugðuÞdu ¼ s;
0
The uniformly distributed delay kernel (27) has been successfully used in a number of different contexts, including models of traffic dynamics with delayed driver response,41 stem cell dynamics,50 time-delayed feedback control,29 and genetic regulation.51 The second example we consider is that of the gamma distribution, which can be written as
The eigenvalues of this general matrix C are given by
N X
ð1
0
for all-to-all coupling with self-feedback. All of the above-considered coupling matrices of unidirectional, bi-directional, and all-to-all couplings with or without self-feedback are special cases of an important class of networks, whose topology is described by the so-called circulant matrices, which are N N matrices of the form
cN
sm hsi ¼
Var ¼
for bi-directional ring coupling with self-feedback, and Gall s;ij ¼
(27)
This distribution has the mean time delay
(22)
for uni-directional ring coupling with self-feedback, 0 1 1 1 0 0 1 B C B1 1 1 0 0C B C B C B C 0 1 1 1 0 1 B C Gbi C; . . . . s ¼ B . 3 B .. .. .. .. .. C B C B C B0 0 0 0 1C @ A 1 0 0 0 1 1 2pk 1 þ 2 cos ; k ¼ 0; …; N 1; k ¼ 3 N
0 c 1 B c2 B B . B C ¼ B .. B B @ cN1
8 < 1 for s q u s þ q; gðuÞ ¼ 2q : 0 elsewhere:
ð1 0
2 ðu sm Þ gðuÞdu ¼
p : a2
(26)
where xk ¼ expð2pik=NÞ. In terms of the distributed delay kernel, we consider two practically important choices of the delay kernel: the uniformly distributed kernel and the gamma distributed delay kernel. The uniformly distributed delay kernel is given by
The gamma distributed delay kernel (30) was originally analysed in models of population dynamics,52–54 and has subsequently been used to study machine tool vibrations,55 intracellular dynamics of HIV infection,56 traffic dynamics with delayed driver response,57 and control of objects over wireless communication networks.58
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 157.211.3.15 On: Mon, 12 Jan 2015 02:35:02
043117-6
€ll Kyrychko, Blyuss, and Scho
Chaos 24, 043117 (2014)
Once the distributed delay kernel is fixed, one can explicitly compute the functions Fc;s and FLc;s as follows. For the uniform distribution (27), they are given by 1 sinðqXÞcosðh XsÞ; qX 1 Fs ðX; hÞ ¼ sinðqXÞsinðh XsÞ; qX
Fc ðX; hÞ ¼
and the case of a discrete time delay can be recovered by setting q ¼ 0. For the weak distribution kernel (30) with p ¼ 1, we have a cos h þ X sin h ; a 2 þ X2 a sin h X cos h ; Fs ðX; hÞ ¼ a a 2 þ X2
Fc ðX; hÞ ¼ a
and similarly, for the strong delay kernel (30) with p ¼ 2 Fc ðX; hÞ ¼ a2 Fs ðX; hÞ ¼ a
ða2 X2 Þcos h þ 2aX sin h
ð a 2 þ X 2 Þ2 2Þ ð 2 2 a X sin h 2aX cos h ð a 2 þ X 2 Þ2
; :
In the same way, Laplace transforms with the modified kernels give FLc ðX; h; KÞ ¼
eKs ½KeKq cos h Xðs qÞ 2 2 2qðK þ X Þ KeKq cos h Xðs þ qÞ þ XeKq sin h Xðs qÞ XeKq sin h Xðs þ qÞ ;
As a first step in the analysis of stability of the fully synchronized solution (12), one should note that this state does not always exist for arbitrary values of the coupling strength and phase and parameters of the delay distribution. The regions of parameter space where this state exists are determined by the condition r02 ¼ k þ Kl½Fc ðX; hÞ cos h 0. Figures 2(a) and 2(b) illustrate such regions for a uniform delay distribution with different choices of the mean time delay, and suggest a natural conclusion that increasing the linear rate k increases the range of values of the coupling strength K for which a fully synchronized solution exists. Once the parameter regions where the fully synchronized state exists have been identified, we choose particular values of the coupling strength K and the mean time delay s, and then compute the master stability function Re ½Kmax ðw; bÞ by solving Eq. (18), as shown in Figures 2(c)–2(f). Stability for individual coupling topologies is verified by checking that Re ½Kmax ðw; bÞ < 0 for all transverse eigenvalues k of the coupling matrix G with Kk ¼ w þ ib. These figures suggest that for the same value of the mean time delay, a uniform distribution with a larger width of the distribution results in the same or a slightly larger region of stable synchronization in the ðw; bÞ plane. Note, however, that this region might also shrink, according to the choice of the mean delay. In the particular case of mean time delay s ¼ 2p, all coupling topologies result in a stable synchronization manifold regardless of the width of the distribution q. On the other hand, for s ¼ 0:52p, the synchronization
Kq eKs Ke sin h Xðs qÞ 2 2Þ ð 2q K þ X KeKq sin h Xðs þ qÞ þ XeKq cos h Xðs þ qÞ XeKq cos h Xðs qÞ ;
FLs ðX; h; KÞ ¼
for the uniform distribution kernel, FLc ðX; h; KÞ ¼ a FLs ðX; h; KÞ ¼ a
ða þ KÞcos h þ X sin h ða þ K Þ2 þ X 2 ða þ KÞsin h X cos h ða þ K Þ2 þ X 2
; ;
for the weak distribution kernel, and h i ða þ KÞ2 X2 cos h þ 2ða þ KÞX sin h ; FLc ðX; h; KÞ ¼ a2 h i2 ða þ KÞ2 þ X2 h i ða þ KÞ2 X2 sin h 2ða þ KÞX cos h ; Fs ðX; hÞ ¼ a2 h i2 ða þ KÞ2 þ X2 for the strong distribution kernel.
FIG. 2. (a) and (b) Existence of synchronized solution and (c)–(f) master stability function Re ½Kmax ðw; bÞ for coupled Stuart-Landau oscillators with uniform delay coupling. Parameter values are x ¼ 1, c ¼ 0, h ¼ 0, and l ¼ 1. (a) s ¼ 2p, (b) s ¼ 0:52p; k ¼ 0:05 (black solid), k ¼ 0:1 (blue dashed), k ¼ 0:2 (red dotted), fully synchronized state exists below the corresponding curves. (c) and (d) Master stability function, k ¼ 0:1, K ¼ 0.3, s ¼ 2p, (c) q ¼ 0, (d) q ¼ 1:49. (e) and (f) Master stability function, k ¼ 0:1, K ¼ 0.08, s ¼ 0:52p, (e) q ¼ 0, (f) q ¼ 0:5p. All eigenvalues of the coupling matrix lie on or inside the black circle.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 157.211.3.15 On: Mon, 12 Jan 2015 02:35:02
043117-7
€ll Kyrychko, Blyuss, and Scho
Chaos 24, 043117 (2014)
studying other types of synchronization, with one particularly important type being that of cluster and splay states. In this case, all oscillators still have the same amplitude r0;m and the same common frequency Xm, but rather than having equal phases, now they have constant differences in their phases rk ðtÞ ¼ r0;m ;
/k ðtÞ ¼ Xm t þ kD/m ;
k ¼ 1; …; N; (32)
with D/m ¼ 2pm=N, where m ¼ 0; …; N 1. The case m ¼ 0 corresponds to identical (or in-phase) synchronization considered earlier, while m ¼ 1; 2; …; N 1 correspond to cluster and splay states. The cluster number Nc ¼ lcm ðm; NÞ=m, where lcmð; Þ denotes the least common multiple, determines how many clusters of oscillators exist, with Nc ¼ N corresponding to a splay state.21,59 As it has been shown in Ref. 22, the cluster solution (32) exists if and only if the following conditions hold N X
Gkj cos ½ðj kÞD/m ¼ Gcm ¼ const;
j¼1 N X
Gkj sin ½ðj kÞD/m ¼ Gsm ¼ const;
(33)
j¼1
FIG. 3. Same as Fig. 2 with gamma distributed delay coupling. Parameter values are x ¼ 1, c ¼ 0, h ¼ 0, and l ¼ 1. (a) weak kernel (p ¼ 1), k ¼ 0:1 (black solid), k ¼ 0.25 (blue dashed), k ¼ 0:4 (red dotted), fully synchronized state exists below the corresponding curve; (b) strong kernel (p ¼ 2), k ¼ 0:15 (black solid), k ¼ 0.25 (blue dashed), k ¼ 0:35 (red dotted), fully synchronized state exists below the corresponding curve. (c) and (d) weak kernel (p ¼ 1), k ¼ 0.25, K ¼ 0.5, (c) a ¼ 0:8. (d) a ¼ 1. (e) and (f) strong kernel (p ¼ 2), (e) a ¼ 1:5, (f) a ¼ 3. All eigenvalues of the coupling matrix lie on or inside the black circle.
independently of k. Although these conditions can be satisfied by different kinds of matrices (see Ref. 22 for a detailed discussion), importantly, they are satisfied by all circulant matrices, of which many standard network topologies are special cases. Provided that the above conditions hold, the common radius r0;m and the common frequency Xm are determined by the following equations: 2 ¼ k þ Kl½Gcm ðXm ; hÞ cos h; r0;m
manifold is always unstable for uni-directional coupling (both with and without self-feedback) and is always stable for bi-directional and all-to-all couplings without and with self-feedback. The values of the coupling strength K were taken to be the same as in the earlier work on synchronization of coupled Stuart-Landau oscillators21,22 to illustrate the role played by the width of the delay distribution. Figure 3 illustrates the computation of regions of existence of the fully synchronized solution, and the master stability function for the gamma delay distribution in the case of p ¼ 1 and p ¼ 2. Similarly to the case of the uniform delay distribution, increasing the bifurcation parameter k leads to an increase in the region of parameters, where the fully synchronized solutions exist. For both weak (p ¼ 1) and strong (p ¼ 2) gamma distribution kernels, decreasing the average time delay (which is equivalent to increasing the parameter a in the gamma distribution) leads to a wider region of stability. Unlike the case of the uniform distribution kernel, varying a can lead to (de)stabilization of certain coupling topologies: while for small values of a only the bidirectional and all-to-all couplings can be stable, as a increases the uni-directional coupling is also stabilized. V. STABILITY OF CLUSTER AND SPLAY STATES
So far, we have considered identical or zero-lag synchronization, where all oscillators have exactly the same amplitudes and phases. It is possible to use a similar approach for
(34)
2 Xm ¼ x cr0;m þ Kl½Gsm ðXm ; hÞ sin h;
where Gcm ðXm ; hÞ ¼ Gcm Fc ðXm ; hÞ Gsm Fs ðXm ; hÞ; Gsm ðXm ; hÞ ¼ Gsm Fc ðXm ; hÞ þ Gcm Fs ðXm ; hÞ: Using the ansatz rk ðtÞ ¼ r0;m ½1 þ drk ðtÞ; /k ðtÞ ¼ Xm t þ kD/m þ d/k ðtÞ, the variational equation for linearization near the cluster state (32) can be found as follows: n_ k ¼ ðJ0;m KlGÞnk þ K Gkj
ð1
N X j¼1
gðt0 ÞRj;m ðt0 Þnj ðt t0 Þdt0 ; k ¼ 1; …; N;
(35)
0
where nk ¼ ðdrk ; d/k ÞT , and we have also introduced the matrices 0 J0;m ¼ @
G¼
0
Rj;m ðt Þ ¼
2 2r0;m
0
2 2cr0;m
0
1 A;
Gcm ðXm ; hÞ Gsm ðXm ; hÞ Gsm ðXm ; hÞ cos Um j sin Um j
Gcm ðXm ; hÞ ! sin Um j ; cos Um j
! ;
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 157.211.3.15 On: Mon, 12 Jan 2015 02:35:02
043117-8
€ll Kyrychko, Blyuss, and Scho
Chaos 24, 043117 (2014)
0 and the phase difference Um j ¼ Xm t þ ðj kÞD/m þ h. Using the same approach as before, one can recast this equation in the form ð1 _n ¼ IN ðJ0;m KlGÞn þ K gðt0 Þ½G Rj;m n ðt t0 Þdt0 ; 0
(36) with n ¼ ðn1 ; …; nN ÞT . To make further progress in the analysis of stability of the cluster state, one has to diagonalize the matrix G and study Floquet exponents corresponding to transverse eigenvalues of this matrix. However, unlike the case of identical synchronization m ¼ 0, this is, in general, impossible, unless the matrix G has certain properties which make the matrix Rj;m independent of k. This is the case for all topologies considered in Sec. IV, which are represented by circulant matrices. To illustrate the computation of the master stability function for cluster states, we consider the network of the StuartLandau oscillators with uni-directional ring coupling: Gk;kþ1 ¼ GN;1 ¼ 1 ¼ l, and all other Gkj ¼ 0. In this case, we 0 have Um j ¼ Xm t þ D/m þ h, and the matrix Rj;m reduces to Rj;m ðt0 Þ ¼ R1 ðt0 Þ ¼
cosðh Xm t0 þ D/m Þ sinðh Xm t0 þ D/m Þ sinðh Xm t0 þ D/m Þ
cosðh Xm t0 þ D/m Þ
! :
Since this matrix does not depend on k, one can use the same approach as in the case of a fully synchronized solution to diagonalize the matrix G, which results in the blockdiagonalized variational equation ð1 f_ k ðtÞ ¼ ðJ0;m KlGÞfk ðtÞ þ Kk gðt0 ÞR1 ðt0 Þfðt t0 Þdt0 ;
In the case of a single discrete time delay gðuÞ ¼ dðu sÞ, this equation reduces to a case investigated by Choe et al.21,22 Figure 4 illustrates the computation of the master stability function of the cluster state for uni-directional ring coupling and uniform or gamma delay distributions. For both of these delay distributions, an increase in the number of oscillators appears to increase the range of coupling strengths, for which the cluster state exists, and also proportionally increase the region of stability of this state. In the case of the gamma distributed delay kernel, one can note that whenever it exists, the cluster state is stable for lower values of the coupling strength across the full range of admissible values of parameter a. VI. DISCUSSION
In this paper, we have shown how the framework of the master stability function for the computation of stability of zero-lag and cluster synchronized states in systems of coupled oscillators can be generalized to the case of distributed delay coupling. To illustrate how the master stability function framework can be used for studying the stability of synchronized states, we have considered the example of a network of coupled Stuart-Landau oscillators, which is a normal form of a supercritical Hopf bifurcation. In this example, computation of the master stability function, i.e., the Floquet exponents, reduces to finding the corresponding eigenvalues of a characteristic equation, which can be done semi-analytically in terms of system parameters,
0
where the eigenvalues k of the matrix G are explicitly given by k ¼ e2ikp=N ; k ¼ 0; 1; …; N 1. Since all terms in the above equation are independent of time, the Floquet exponents can be found as the eigenvalues K of the characteristic equation ð1 0 0 Kt0 0 gðt ÞR1 ðt Þe dt ¼ 0; det J0;m KlG KI2 þ Kk 0
(37) where I2 is a 2 2 identity matrix. More explicitly, this transcendental equation has the form 2 þ KðlGcm ðXm ; hÞ k FLc ðXm ; h þ D/m ; KÞÞK K2 þ 2½r0;m 2 K½lGcm ðXm ; hÞ k FLc ðXm ; h þ D/m ; KÞ þ 2r0;m
þ cðlGsm ðXm ; hÞ k FLs ðXm ; h þ D/m ; KÞÞ þ K 2 ½lGcm ðXm ; hÞ k FLc ðXm ; h þ D/m ; KÞ2 þ K 2 ½lGsm ðXm ; hÞ k FLs ðXm ; h þ D/m ; KÞ2 þ K 2 ½lGsm ðXm ; hÞ k FLs ðXm ; h þ D/m ; KÞ2 ¼ 0: (38)
FIG. 4. Master stability function for a splay state with m ¼ 1 for a ring of uni-directionally coupled Stuart-Landau oscillators. Parameter values are x ¼ 1, c ¼ 0, h ¼ 0, and l ¼ 1. (a) and (b) uniform delay distribution, s ¼ 2p; k ¼ 0:1, (a) N ¼ 4, (b) N ¼ 10. (c) and (d) weak kernel (p ¼ 1), k ¼ 0.25, (c) N ¼ 4, (d) N ¼ 10. (e) and (f) strong kernel (p ¼ 2), k ¼ 0.25, (e) N ¼ 4, (f) N ¼ 10. The splay state does not exist in the white region.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 157.211.3.15 On: Mon, 12 Jan 2015 02:35:02
043117-9
€ll Kyrychko, Blyuss, and Scho
topological eigenvalues of the coupling matrix, and the coupling parameters, where the Laplace transform of the delay kernel enters. Calculations for the case of uniform and gamma distributed delay kernels show that increasing the width of the delay distribution (for a uniform kernel) or reducing the mean time delay (for a gamma distributed delay kernel) leads often to a larger region of stability, thus allowing the fully synchronized state to be stable for a larger class of coupling topologies. Although the analysis presented in this paper has focused upon the case of linear coupling between different oscillators in the network, it is straightforward to generalize our results on the stability of zero-lag and cluster synchronization to systems with nonlinear coupling,60 since the master stability equation is still a linear variational equation involving the appropriate delay distribution kernel. Recent studies of amplitude death in systems of delay-coupled oscillators47,48,61 have highlighted the important role played by the coupling phase in determining the regions of possible suppression of oscillations. Due to substantial analytical headway provided by considering Stuart-Landau oscillators, it would be insightful to use the master stability framework developed in this paper to investigate the effects of the coupling phase on existence and stability of synchronized solutions in networks of such oscillators. ACKNOWLEDGMENTS
This work was partially supported by DFG in the framework of SFB 910: Control of self-organizing nonlinear systems: Theoretical methods and concepts of application. Y.N.K. and K.B.B. gratefully acknowledge the hospitality of the Institut f€ ur Theoretische Physik, TU Berlin, where part of this work was completed. E.S. is grateful to Carolin Wille for helpful discussions. The authors thank the anonymous referees for helpful comments and suggestions that have helped to improve the presentation.
1
R. Albert and A.-L. Barabasi, Rev. Mod. Phys. 74, 47 (2002). A.-L. Barabasi and E. Bonabeau, Sci. Am. 288, 60 (2003). 3 M. E. J. Newman, SIAM Rev. 45, 167 (2003). 4 D. J. Watts and S. H. Strogatz, Nature 393, 440 (1998). 5 S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D. U. Hwang, Phys. Rep. 424, 175 (2006). 6 A. Pikovsky, M. Rosenblum, and J. Kurths, Synchronization: A Universal Concept in Nonlinear Sciences (Cambridge University Press, Cambridge, 2001). 7 L. M. Pecora and T. L. Carroll, Phys. Rev. Lett. 80, 2109 (1998). 8 M. G. Earl and S. H. Strogatz, Phys. Rev. E 67, 036204 (2003). 9 O. D’Huys, R. Vicente, T. Erneux, J. Danckaert, and I. Fischer, Chaos 18, 037116 (2008). 10 W. Kinzel, A. Englert, G. Reents, M. Zigzag, and I. Kanter, Phys. Rev. E 79, 056207 (2009). 11 O. D’Huys, R. Vicente, J. Danckaert, and I. Fischer, Chaos 20, 043127 (2010). 12 A. Englert, W. Kinzel, Y. Aviad, M. Butkovski, I. Reidler, M. Zigzag, I. Kanter, and M. Rosenbluh, Phys. Rev. Lett. 104, 114102 (2010). 13 V. Flunkert, S. Yanchuk, T. Dahms, and E. Sch€ oll, Phys. Rev. Lett. 105, 254101 (2010). 14 Complex Time-Delay Systems, edited by F. M. Atay (Springer, New York, 2010). 15 A. Englert, S. Heiligenthal, W. Kinzel, and I. Kanter, Phys. Rev. E 83, 046222 (2011). 2
Chaos 24, 043117 (2014) 16
K. Hicke, O. D’Huys, V. Flunkert, E. Sch€ oll, I. Danckaert, and I. Fischer, Phys. Rev. E 83, 056211 (2011). J. Lehnert, T. Dahms, P. H€ ovel, and E. Sch€ oll, Europhys. Lett. 96, 60013 (2011). 18 V. Flunkert, I. Fischer, and E. Sch€ oll, “Dynamics, control and information in delay-coupled systems,” Philos. Trans. R. Soc. London, Ser. A 371, 20120465 (2013). 19 M. C. Soriano, J. Garcıa-Ojalvo, C. R. Mirasso, and I. Fischer, Rev. Mod. Phys. 85, 421 (2013). 20 A. Zakharova, I. Schneider, Y. N. Kyrychko, K. B. Blyuss, A. Koseska, B. Fiedler, and E. Sch€ oll, Europhys. Lett. 104, 50004 (2013). 21 C.-U. Choe, T. Dahms, P. H€ ovel, and E. Sch€ oll, Phys. Rev. E 81, 025205(R) (2010). 22 C.-U. Choe, T. Dahms, P. H€ ovel, and E. Sch€ oll, in Proceedings of the Eighth AIMS International Conference on Dynamical Systems, Differential Equations and Applications (American Institute of Mathematical Sciences, Springfield, 2011), pp. 292–301. 23 C.-U. Choe, H. Jang, V. Flunkert, T. Dahms, P. H€ ovel, and E. Sch€ oll, Dyn. Syst. 28, 15 (2013). 24 T. Dahms, J. Lehnert, and E. Sch€ oll, Phys. Rev. E 86, 016202 (2012). 25 M. Dhamala, V. K. Jirsa, and M. Ding, Phys. Rev. Lett. 92, 074104 (2004). 26 C. R. S. Williams, T. E. Murphy, R. Roy, F. Sorrentino, T. Dahms, and E. Sch€ oll, Phys. Rev. Lett. 110, 064104 (2013). 27 D. Hunt, G. Korniss, and B. K. Szymanski, Phys. Rev. Lett. 105, 068701 (2010). 28 D. Hunt, B. K. Szymanski, and G. Korniss, Phys. Rev. E 86, 056114 (2012). 29 A. Gjurchinovski and V. Urumov, Europhys. Lett. 84, 40013 (2008). 30 A. Gjurchinovski and V. Urumov, Phys. Rev. E 81, 016209 (2010). 31 A. Gjurchinovski, A. Zakharova, and E. Sch€ oll, Phys. Rev. E 89, 032915 (2014). 32 F. Atay, Phys. Rev. Lett. 91, 094101 (2003). 33 S. A. Campbell and R. Jessop, Math. Model. Nat. Phenom. 4, 1 (2009). 34 C. U. Choe, R.-S. Kim, H. Jang, P. H€ ovel, and E. Sch€ oll, Int. J. Dyn. Control 2, 2 (2014). 35 G. Kiss and B. Krauskopf, Dyn. Syst. 26, 85 (2011). 36 A. Gjurchinovski, T. J€ ungling, V. Urumov, and E. Sch€ oll, Phys. Rev. E 88, 032912 (2013). 37 W. Michiels, V. Van Assche, and S.-I. Niculescu, IEEE Trans. Autom. Control 50, 493 (2005). 38 K. B. Blyuss and Y. N. Kyrychko, Bull. Math. Biol. 72, 490 (2010). 39 T. Faria and S. Trofimchuk, Nonlinearity 23, 2457 (2010). 40 S. A. Gourley and J. W.-H. So, Proc. R. Soc. Edinburgh 133, 527 (2003). 41 R. Sipahi, F. M. Atay, and S.-I. Niculescu, SIAM J. Appl. Math. 68, 738 (2007). 42 C. W. Eurich, A. Thiel, and L. Fahse, Phys. Rev. Lett. 94, 158104 (2005). 43 R. Vicente, L. L. Gollo, C. R. Mirasso, I. Fischer, and P. Gordon, Proc. Natl. Acad. Sci. U.S.A. 105, 17157 (2008). 44 C. Wille, J. Lehnert, and E. Sch€ oll, Phys. Rev. E 90, 032908 (2014). 45 I.-C. Morarescu, W. Michiels, and M. Jungers, in Proceedings of the American Control Conference, Washington (2013). 46 B. Fiedler, V. Flunkert, M. Georgi, P. H€ ovel, and E. Sch€ oll, Phys. Rev. Lett. 98, 114101 (2007). 47 Y. N. Kyrychko, K. B. Blyuss, and E. Sch€ oll, Eur. Phys. J. B 84, 307 (2011). 48 Y. N. Kyrychko, K. B. Blyuss, and E. Sch€ oll, Philos. Trans. R. Soc. London, Ser. A 371, 20120466 (2013). 49 K. Pyragas and T. Pyragiene, Phys. Rev. E 78, 046217 (2008). 50 M. Adimy, F. Crauste, and S. Ruan, Nonlinear Anal.: Real World Appl. 6, 651 (2005). 51 M. Barrio, K. Burrage, A. Leier, and T. Tian, PLoS Comput. Biol. 2, e117 (2006). 52 S. P. Blythe, R. M. Nisbet, W. S. C. Gurney, and N. MacDonald, J. Math. Anal. Appl. 109, 388 (1985). 53 K. Cooke and Z. Grossman, J. Math. Anal. Appl. 86, 592 (1982). 54 J. M. Cushing, in Mathematics of Biology, edited by M. Iannelli (Springer, New York, 2011). 55 G. Stepan, in Dynamics and Chaos in Manufacturing Process, edited by F. C. Moon (Wiley, New York, 1998). 56 J. E. Mittler, B. Sulzer, A. U. Neumann, and A. S. Perelson, Math. Biosci. 152, 143 (1998). 57 R. Sipahi, S.-I. Niculescu, and F. M. Atay, in Proceedings of the American Control Conference, New York (2007). 17
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 157.211.3.15 On: Mon, 12 Jan 2015 02:35:02
043117-10 58
€ ll Kyrychko, Blyuss, and Scho
O. Roesch, H. Roth, and S.-I. Niculescu, in Proceedings of the IEEE International Conference on Mechatronics and Automation, Niagara Falls (2005). 59 R. Zillmer, R. Livi, A. Politi, and A. Torcini, Phys. Rev. E 76, 046102 (2007).
Chaos 24, 043117 (2014) 60
61
O. V. Popovych, C. Hauptmann, and P. A. Tass, Phys. Rev. Lett. 94, 164102 (2005). W. Zou, J. Lu, Y. Tang, C. Zhang, and J. Kurths, Phys. Rev. E 84, 066208 (2011).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 157.211.3.15 On: Mon, 12 Jan 2015 02:35:02
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 157.211.3.15 On: Mon, 12 Jan 2015 02:35:02
CHAOS 24, 043117 (2014)
Synchronization of networks of oscillators with distributed delay coupling €ll2 Y. N. Kyrychko,1,a) K. B. Blyuss,1 and E. Scho 1
Department of Mathematics, University of Sussex, Falmer, Brighton BN1 9QH, United Kingdom Institut f€ ur Theoretische Physik, Technische Universit€ at Berlin, Hardenbergstr. 36, 10623 Berlin, Germany
2
(Received 16 June 2014; accepted 9 October 2014; published online 24 October 2014) This paper studies the stability of synchronized states in networks, where couplings between nodes are characterized by some distributed time delay, and develops a generalized master stability function approach. Using a generic example of Stuart-Landau oscillators, it is shown how the stability of synchronized solutions in networks with distributed delay coupling can be determined through a semi-analytic computation of Floquet exponents. The analysis of stability of fully synchronized and of cluster or splay states is illustrated for several practically important choices of C 2014 AIP Publishing LLC. delay distributions and network topologies. V [http://dx.doi.org/10.1063/1.4898771] Many coupled real-life systems exhibit complex dynamical regimes, including the emergence of coherent dynamics known as synchronization. An important practical issue is understanding for what types of networks the synchronization can occur and is stable, which determines whether or not it can be observed experimentally. Quite often, connections between different elements of the network are non-instantaneous, and one has to explicitly include time delays in these connections in the analysis of synchronization. In this paper, we show how one can analyze stability of different types of synchronization in networks with distributed delay coupling, which represents a realistic situation in which the time delays themselves may not be constant or fixed. Our results suggest that the parameter space and network topologies for which a stable synchronization occurs are affected not only by the mean time delay but also by the width of the delay distribution.
I. INTRODUCTION
In the last decade, there has been a substantial growth of research interest in the dynamics of coupled systems, ranging from a few elements to large networks.1–4 From a perspective of potential applications, one of the central research questions for such systems is the emergence and stability of different types of collective dynamics and, in particular, synchronization.5,6 In order to study the stability of synchronization, Pecora and Carroll7 put forward a master stability function (MSF) approach that allows one to separate the local dynamics of individual nodes from the network topology, which is achieved by a proper diagonalization of the matrix representing the full network dynamics. An important aspect of network dynamics concerns the fact that interactions between nodes are often noninstantaneous due to a finite speed of signal propagation. In order to account for this feature, one has to explicitly include time delays in the analysis, and this is known to have a a)
Electronic mail: [email protected]
1054-1500/2014/24(4)/043117/10/$30.00
significant impact on network behaviour.8–20 Some work has been done on the extension of the master stability function to networks with time-delayed coupling10,21–26 but so far it has only included single constant time delays. Hunt et al.27,28 have considered complex networks of linearly coupled nodes with time delays and noise and established conditions for synchronizability of such networks. At the same time, in many realistic systems, time delays may be not constant and may either vary depending on the values of system variables, or just not be explicitly known,29–31 and in these cases, the standard methodology of considering one or several constant time delays is not sufficient. To describe such situations mathematically, one can use the formalism of distributed time delays, where the time delay is represented through an integral kernel describing a particular delay distribution.32–35 Distributed time delays have been successfully used to describe situations when only an approximate value of time delay is known in engineering experiments,35–37 for modelling distributions of waiting times in epidemiological models,38 maturation periods in population and ecological models,39,40 as well as in models of traffic dynamics,41 and neural systems.42–44 In a recent paper, Morarescu et al.45 have looked at systems with gamma-distributed delay coupling and analyzed the stability of the synchronized equilibria (steady states), which are not affected by the coupling. In this paper, we consider networks of coupled identical dynamical systems with linear distributed delay coupling as represented by some integral kernel. Linearization near the synchronization manifold yields a variational equation, which after block-diagonalization of the coupling matrix reduces to a single complex-valued integro-differential equation, whose maximum Floquet or Lyapunov exponent gives the master stability function in terms of the system parameters. Using the example of a network of Stuart-Landau oscillators, we make further analytical progress by showing how the computation of the master stability function can be reduced to finding solutions of a single transcendental characteristic equation. For particular choices of the delay distribution kernel, all terms in this equation can be found in a closed analytical form in terms of coupling and system parameters. The methodology we develop is quite generic and
24, 043117-1
C 2014 AIP Publishing LLC V
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 157.211.3.15 On: Mon, 12 Jan 2015 02:35:02
043117-2
€ll Kyrychko, Blyuss, and Scho
Chaos 24, 043117 (2014)
can be applied to analyze the stability of fully synchronized and of cluster states in any systems with distributed delay coupling. The outline of this paper is as follows. In Sec. II, we develop a general master stability function formalism for networks of identical coupled elements with distributed delay coupling. Section III illustrates how this general framework can be made more explicit by means of an amplitude-phase representation for the example of coupled Stuart-Landau oscillators. Section IV contains the results of analytical and numerical computations of the master stability function for coupled Stuart-Landau oscillators with uniform and gamma distributed delays and different coupling topologies. In Sec. V, the methodology of the master stability function for systems with distributed delay coupling is extended to the analysis of cluster states representing generalized synchronization. The paper concludes with discussion of results in Sec. VI. II. MASTER STABILITY FUNCTION
We consider a network of N identical dynamical systems with distributed delay coupling x_ k ¼ Fk ½xk ðtÞ þ r
N X
ð1 Gkj Hkj
0
0
0
gðt Þxj ðt t Þdt xk ðtÞ ;
are identical for different nodes, i.e., Hkj ¼ H for k; j ¼ 1; …N. The last two assumptions ensure that all nodes received the same input. With those assumptions, the model (1) simplifies to ð1 N X 0 0 0 Gkj H gðt Þxj ðt t Þdt xk ðtÞ ; x_ k ¼ F½xk ðtÞ þ r k ¼ 1; …; N:
k ¼ 1; …; N;
(1)
(3)
Let us now consider the synchronization manifold, defined as xk ðtÞ xs ðtÞ for all k ¼ 1; …N. The dynamics on the synchronization manifold is given by ð1 0 0 0 gðt Þxs ðt t Þdt xs ðtÞ ; (4) x_ s ¼ F½xs ðtÞ þ rlH 0
where l is the row sum, Eq. (2). It is worth noting that the dynamics of synchronization manifold is independent of the coupling strength in the case of zero row sum l ¼ 0 or instantaneous coupling gðuÞ ¼ dðuÞ, and in all other cases it will also depend on the coupling strength r. To study the stability of the synchronization manifold, we linearize the full system (3) near xk ðtÞ ¼ xs ðtÞ, which yields with xk ðtÞ ¼ xs ðtÞ þ nk ðtÞ,
0
j¼1
0
j¼1
n_ k ðtÞ ¼ J0 ðtÞnk ðtÞ þ r
N X j¼1
m
where xk 2 R is the state vector of each system, Fk ½ðxk ðtÞ describes the local dynamics of the node k, r 2 C is the coupling strength, G is an N N matrix describing the coupling topology and the strength of each link in the network, Hkj is an m m matrix representing the coupling scheme, and gðÞ is a delay distribution kernel, satisfying ð1 gðuÞdu ¼ 1: gðuÞ 0; 0
For gðuÞ ¼ dðuÞ, one recovers instantaneous coupling ðxj xk Þ; for gðuÞ ¼ dðu sÞ, the coupling takes the form of a discrete time delay: ½xj ðt sÞ xk ðtÞ. The matrix G, whose non-zero entries Gij define a relative strength of a link between nodes j and i, is related to the adjacency matrix of the network topology, whose elements are zero or one depending on the existence of the respective link. Next, we discuss a number of assumptions regarding the individual dynamical systems and their interactions in the network. All nodes are assumed to be identical, so that their dynamics is described by a single function Fk ½ ¼ F½. The time delays in the propagation of signals from different nodes are also taken to be the same and described by the delay distribution kernel g(u). Since we are interested in synchronization and its stability, we will further assume that the matrix G has a constant row sum N X
Gkj ¼ l:
(2)
ð1 Gkj H 0
k ¼ 1; …; N;
(5)
where J0 ðtÞ ¼ DF½xs ðtÞ is the Jacobian of F. Introducing a vector n ¼ ðn1 ; n2 ; …; nN ÞT , the above equation can be rewritten as ð1 _n ¼ IN ½J0 ðtÞ rlH n þ r½G H gðt0 Þnðt t0 Þdt0 ; 0
(6) where denotes the Kronecker product, and IN is the N N unity matrix. We will assume that matrix G is diagonalizable, i.e., there exists a unitary transformation U such that UGU1 ¼ diagðl; 1 ; 2 ; …; N1 Þ: Here, the first eigenvalue of the matrix G, which is equal to the row sum 0 l, corresponds to the eigenvector ð1; 1; …; 1Þ describing the dynamics on the synchronization manifold. Hence, this eigenvalue is called the longitudinal eigenvalue of G, whereas all other eigenvalues are called transverse eigenvalues as they correspond to directions transverse to the synchronization manifold. To make further analytical progress, it is instructive to diagonalize the coupling matrix G, which results in a blockdiagonalized variational equation
j¼1
This condition is necessary for the existence of the synchronized solution. Finally, we assume that the coupling schemes
gðt0 Þnj ðt t0 Þdt0 nk ðtÞ ;
f_ k ðtÞ ¼ ½J0 ðtÞ rlHfk ðtÞ þ rk H
ð1
gðt0 Þfk ðt t0 Þdt0 ;
0
k ¼ 0; …; N 1;
(7)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 157.211.3.15 On: Mon, 12 Jan 2015 02:35:02
043117-3
€ll Kyrychko, Blyuss, and Scho
Chaos 24, 043117 (2014)
where k is the k-th eigenvalue of the coupling matrix G. It is worth noting that in this system of equations, the Jacobian J0 ðtÞ and H are the same for all k, and one can consider the variational equation in dependence on the complex parameter w þ ib ð1 _fðtÞ ¼ ½J0 ðtÞ rlHfðtÞ þ ðw þ ibÞH gðt0 Þfðt t0 Þdt0 : 0
(8) Computation of the maximum Lyapunov exponent (or real part of the Floquet exponent) for this equation yields the master stability function
Introducing amplitude and phase variables as rk ¼ jzk j and /k ¼ argðzk Þ, Eq. (9) can be recast as a system of 2N real equations ð1 N X gðt0 Þrj ðt t0 Þcos r_ k ¼ k rk2 ðtÞ rk ðtÞ þ K Gkj 0
j¼1
i /j ðt t0 Þ /k ðtÞ þ h dt0 rk ðtÞcos h ; h
ð1 N X rj ðt t 0 Þ sin Gkj gðt0 Þ rk ðtÞ 0 j¼1 0 0 ½/j ðt t Þ /k ðtÞ þ hdt sin h ; k ¼ 1; …; N:
/_ k ¼ x crk2 þ K
Re½Kmax ðw; b; rÞ
(10)
and the analysis of stability of the synchronization manifold reduces to checking that for all transverse eigenvalues of the coupling matrix rk ¼ w þ ib; k ¼ 1; …; N 1, one has Re[Kmax ðw; b; rÞ < 0. A positive value of the Lyapunov exponent (real part of the Floquet exponent) associated with the longitudinal eigenvalue 0 indicates that the synchronized solution is chaotic (unstable). The significant advantage of this approach lies in the separation of the actual dynamics from the topology of the coupling, as the computation of Kmax ðw; b; rÞ can be done once for Eq. (8) independently of the coupling matrix. It is worth noting that although the computation of the master stability function can be performed independently of the coupling topology, it has to be done separately for each value of the coupling strength r. The reason for this is the fact that the row sum l is, in general, non-zero, and also the synchronized state xs ðtÞ itself depends on the coupling strength due to the delayed nature of the coupling, as follows from Eq. (4). The above approach can be generalized in a straightforward manner to more complex forms of the coupling function H. III. COUPLED STUART-LANDAU OSCILLATORS
In order to illustrate the analysis of the stability of synchronization, we consider an array of N identical diffusively coupled Stuart-Landau oscillators (k ¼ 1; …; N) z_ k ðtÞ ¼ ðk þ ixÞzk ðtÞ ð1 þ icÞjzk ðtÞj2 zk ðtÞ ð1 N X 0 0 0 þr Gkj gðt Þzj ðt t Þdt zk ðtÞ ; j¼1
(9)
0
K 2 Rþ ;
rk ðtÞ ¼ r0 ;
/k ¼ Xt;
k ¼ 1; …; N;
(11)
with the common radius r0 and the common frequency X of oscillations being determined by the solutions of the following equations: r02 ¼ k þ Kl½Fc ðX; hÞ cos h;
(12)
X ¼ x cr02 þ Kl½Fs ðX; hÞ sin h; where we have introduced the auxiliary quantities ð1 gðt0 Þ cosðat0 þ bÞdt0 ; Fc ða; bÞ ¼ 0 ð1 gðt0 Þ sinðat0 þ bÞdt0 : Fs ða; bÞ ¼
(13)
0
To analyse the stability of the fully synchronized solution (11), we use a slightly modified ansatz for small perturbations around this solution, namely, rk ðtÞ ¼ r0 ½1 þ drk ðtÞ;
/k ðtÞ ¼ Xt þ d/k ðtÞ;
k ¼ 1; …; N; which yields a variational equation of the form n_ k ¼ ðJ0 KlR0 Þnk þ K
N X
ð1 Gkj
k ¼ 1; …; N;
(14)
where nk ¼ ðdrk ; d/k ÞT , and the matrices J0 ; R0 and R1 are given by
h 2 R;
where K and h are the strength and the phase of coupling, respectively. Such complex-valued couplings, which are equivalent to a rotational matrix H if zk is written in terms of real and imaginary part, have been shown to be important in overcoming the odd-number limitation of time-delayed feedback control,46 in controlling amplitude death,47,48 and in anticipating chaos synchronization.49
gðt0 ÞR1 ðt0 Þnj ðt t0 Þdt0 ;
0
j¼1
which is a prototype of dynamics near a supercritical Hopf bifurcation. Here, zk 2 C, k, x 6¼ 0, and c are real constants, where x is the intrinsic oscillation frequency, and r ¼ Keih ;
The fully synchronized solution of the system (10) is given by
J0 ¼ R0 ¼ R1 ðt0 Þ ¼
2r02
0
2cr02
0
! ;
Fc ðX; hÞ Fs ðX; hÞ
; Fs ðX; hÞ Fc ðX; hÞ cosðh Xt0 Þ sinðh Xt0 Þ sinðh Xt0 Þ
cosðh Xt0 Þ
:
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 157.211.3.15 On: Mon, 12 Jan 2015 02:35:02
043117-4
€ll Kyrychko, Blyuss, and Scho
Chaos 24, 043117 (2014)
Fc ða; bÞ ¼ FLc ða; b; 0Þ;
Equivalently, the above system can be rewritten as n_ ¼ IN ðJ0 KlR0 Þn ð1 gðt0 Þ½G R1 ðt0 Þ n ðt t0 Þdt0 ; þK
(15)
Fs ða; bÞ ¼ FLs ða; b; 0Þ:
Equation (18) generalises an earlier result of Choe et al.,21,22 which considered the master stability function for the case of a single discrete time delay.
0
with the 2N-dimensional vector n ¼ ðn1 ; …; nN ÞT . Diagonalizing the matrix G results in a block-diagonalized variational equation ð1 gðt0 ÞR1 ðt0 Þfk ðt t0 Þdt0 : f_ k ðtÞ ¼ ðJ0 KlR0 Þfk ðtÞ þ Kk 0
(16) Since in these equations, none of the matrices depends explicitly on time t, it is possible to find Floquet exponents as the eigenvalues K 2 C of the following characteristic equation ð1 gðt0 ÞR1 ðt0 ÞeKt0 dt0 ¼ 0; det J0 KlR0 KI2 þ Kk 0
(17) where I2 is a 2 2 identity matrix. More explicitly, this transcendental equation has the form K
2
þ 2½r02 þ KðlFc ðX; hÞ k FLc ðX; h; KÞÞK þ 2r02 K½lFc ðX; hÞ k FLc ðX; h; KÞ þ2cr02 K½lFs ðX; hÞ k FLs ðX; h; KÞ þ K 2 ½lFc ðX; hÞ k FLc ðX; h; KÞ2 þK 2 ½lFs ðX; hÞ k FLs ðX; h; KÞ2 ¼ 0;
(18)
where by analogy with (13), we have introduced the quantities ð1 0 gðt0 Þ cosðat0 þ bÞezt dt0 ; FLc ða; b; zÞ ¼ 0 ð1 0 gðt0 Þ sinðat0 þ bÞezt dt0 ; FLs ða; b; zÞ ¼
IV. COMPUTATION OF MASTER STABILITY FUNCTION
In order to illustrate how the master stability function can be used for the analysis of stability of a synchronization manifold, we first have to specify a particular topology of the network G, as well as the distributed delay kernel g(u). We will consider a number of different coupling topologies, each of which is characterized by a particular spectrum of eigenvalues fk g, as shown in Fig. 1. Assuming that the constant row sum l is different from zero, we normalize matrices G by dividing all elements by l, which results in a row sum equal to unity. The coupling matrices G and their eigenvalues are then given by 0 1 0 1 0 0 0 B0 0 1 0 0C B C B. . . . .. C uni 2pik=N B . . . . ; G ¼ B . . . . . C C; k ¼ e B C @0 0 0 0 1A 1
0 0
0 (19)
for uni-directional ring coupling, 0 1 0 1 0 0 1 B C B1 0 1 0 0C B C B0 1 0 1 0C B C 1 Gbi ¼ B .. C .. .. .. .. C; 2B B . . . . . C B C B0 0 0 0 1C @ A 1 0 0 0 0
k ¼ cos
2pk ; N
k ¼ 0; …; N 1;
(20)
for bi-directional ring coupling, and
0
where the superscript L refers to these integrals representing the Laplace transform of modified kernels gðsÞ cosðas þ bÞ and gðsÞ sinðas þ bÞ. Comparing these expressions to (13) yields the following relations:
0
k ¼ 0; …; N 1;
Gall ij
8
1 (known as strong delay kernel in the case p ¼ 2), the biggest influence on the coupling at any moment of time t is from the values of zk at t ðp 1Þ=a. The delay distribution (30) has the mean time delay ð1 p (31) ugðuÞdu ¼ ; sm ¼ a 0 and the variance
ð1mÞþN mod N cm x k ;
Var ¼
m¼1
k ¼ 0; …; N 1;
ð1
gðuÞ ¼
k ¼ c1 þ cN xk þ cN1 x2k þ þ c2 xN1 k ¼
ugðuÞdu ¼ s;
0
The uniformly distributed delay kernel (27) has been successfully used in a number of different contexts, including models of traffic dynamics with delayed driver response,41 stem cell dynamics,50 time-delayed feedback control,29 and genetic regulation.51 The second example we consider is that of the gamma distribution, which can be written as
The eigenvalues of this general matrix C are given by
N X
ð1
0
for all-to-all coupling with self-feedback. All of the above-considered coupling matrices of unidirectional, bi-directional, and all-to-all couplings with or without self-feedback are special cases of an important class of networks, whose topology is described by the so-called circulant matrices, which are N N matrices of the form
cN
sm hsi ¼
Var ¼
for bi-directional ring coupling with self-feedback, and Gall s;ij ¼
(27)
This distribution has the mean time delay
(22)
for uni-directional ring coupling with self-feedback, 0 1 1 1 0 0 1 B C B1 1 1 0 0C B C B C B C 0 1 1 1 0 1 B C Gbi C; . . . . s ¼ B . 3 B .. .. .. .. .. C B C B C B0 0 0 0 1C @ A 1 0 0 0 1 1 2pk 1 þ 2 cos ; k ¼ 0; …; N 1; k ¼ 3 N
0 c 1 B c2 B B . B C ¼ B .. B B @ cN1
8 < 1 for s q u s þ q; gðuÞ ¼ 2q : 0 elsewhere:
ð1 0
2 ðu sm Þ gðuÞdu ¼
p : a2
(26)
where xk ¼ expð2pik=NÞ. In terms of the distributed delay kernel, we consider two practically important choices of the delay kernel: the uniformly distributed kernel and the gamma distributed delay kernel. The uniformly distributed delay kernel is given by
The gamma distributed delay kernel (30) was originally analysed in models of population dynamics,52–54 and has subsequently been used to study machine tool vibrations,55 intracellular dynamics of HIV infection,56 traffic dynamics with delayed driver response,57 and control of objects over wireless communication networks.58
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 157.211.3.15 On: Mon, 12 Jan 2015 02:35:02
043117-6
€ll Kyrychko, Blyuss, and Scho
Chaos 24, 043117 (2014)
Once the distributed delay kernel is fixed, one can explicitly compute the functions Fc;s and FLc;s as follows. For the uniform distribution (27), they are given by 1 sinðqXÞcosðh XsÞ; qX 1 Fs ðX; hÞ ¼ sinðqXÞsinðh XsÞ; qX
Fc ðX; hÞ ¼
and the case of a discrete time delay can be recovered by setting q ¼ 0. For the weak distribution kernel (30) with p ¼ 1, we have a cos h þ X sin h ; a 2 þ X2 a sin h X cos h ; Fs ðX; hÞ ¼ a a 2 þ X2
Fc ðX; hÞ ¼ a
and similarly, for the strong delay kernel (30) with p ¼ 2 Fc ðX; hÞ ¼ a2 Fs ðX; hÞ ¼ a
ða2 X2 Þcos h þ 2aX sin h
ð a 2 þ X 2 Þ2 2Þ ð 2 2 a X sin h 2aX cos h ð a 2 þ X 2 Þ2
; :
In the same way, Laplace transforms with the modified kernels give FLc ðX; h; KÞ ¼
eKs ½KeKq cos h Xðs qÞ 2 2 2qðK þ X Þ KeKq cos h Xðs þ qÞ þ XeKq sin h Xðs qÞ XeKq sin h Xðs þ qÞ ;
As a first step in the analysis of stability of the fully synchronized solution (12), one should note that this state does not always exist for arbitrary values of the coupling strength and phase and parameters of the delay distribution. The regions of parameter space where this state exists are determined by the condition r02 ¼ k þ Kl½Fc ðX; hÞ cos h 0. Figures 2(a) and 2(b) illustrate such regions for a uniform delay distribution with different choices of the mean time delay, and suggest a natural conclusion that increasing the linear rate k increases the range of values of the coupling strength K for which a fully synchronized solution exists. Once the parameter regions where the fully synchronized state exists have been identified, we choose particular values of the coupling strength K and the mean time delay s, and then compute the master stability function Re ½Kmax ðw; bÞ by solving Eq. (18), as shown in Figures 2(c)–2(f). Stability for individual coupling topologies is verified by checking that Re ½Kmax ðw; bÞ < 0 for all transverse eigenvalues k of the coupling matrix G with Kk ¼ w þ ib. These figures suggest that for the same value of the mean time delay, a uniform distribution with a larger width of the distribution results in the same or a slightly larger region of stable synchronization in the ðw; bÞ plane. Note, however, that this region might also shrink, according to the choice of the mean delay. In the particular case of mean time delay s ¼ 2p, all coupling topologies result in a stable synchronization manifold regardless of the width of the distribution q. On the other hand, for s ¼ 0:52p, the synchronization
Kq eKs Ke sin h Xðs qÞ 2 2Þ ð 2q K þ X KeKq sin h Xðs þ qÞ þ XeKq cos h Xðs þ qÞ XeKq cos h Xðs qÞ ;
FLs ðX; h; KÞ ¼
for the uniform distribution kernel, FLc ðX; h; KÞ ¼ a FLs ðX; h; KÞ ¼ a
ða þ KÞcos h þ X sin h ða þ K Þ2 þ X 2 ða þ KÞsin h X cos h ða þ K Þ2 þ X 2
; ;
for the weak distribution kernel, and h i ða þ KÞ2 X2 cos h þ 2ða þ KÞX sin h ; FLc ðX; h; KÞ ¼ a2 h i2 ða þ KÞ2 þ X2 h i ða þ KÞ2 X2 sin h 2ða þ KÞX cos h ; Fs ðX; hÞ ¼ a2 h i2 ða þ KÞ2 þ X2 for the strong distribution kernel.
FIG. 2. (a) and (b) Existence of synchronized solution and (c)–(f) master stability function Re ½Kmax ðw; bÞ for coupled Stuart-Landau oscillators with uniform delay coupling. Parameter values are x ¼ 1, c ¼ 0, h ¼ 0, and l ¼ 1. (a) s ¼ 2p, (b) s ¼ 0:52p; k ¼ 0:05 (black solid), k ¼ 0:1 (blue dashed), k ¼ 0:2 (red dotted), fully synchronized state exists below the corresponding curves. (c) and (d) Master stability function, k ¼ 0:1, K ¼ 0.3, s ¼ 2p, (c) q ¼ 0, (d) q ¼ 1:49. (e) and (f) Master stability function, k ¼ 0:1, K ¼ 0.08, s ¼ 0:52p, (e) q ¼ 0, (f) q ¼ 0:5p. All eigenvalues of the coupling matrix lie on or inside the black circle.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 157.211.3.15 On: Mon, 12 Jan 2015 02:35:02
043117-7
€ll Kyrychko, Blyuss, and Scho
Chaos 24, 043117 (2014)
studying other types of synchronization, with one particularly important type being that of cluster and splay states. In this case, all oscillators still have the same amplitude r0;m and the same common frequency Xm, but rather than having equal phases, now they have constant differences in their phases rk ðtÞ ¼ r0;m ;
/k ðtÞ ¼ Xm t þ kD/m ;
k ¼ 1; …; N; (32)
with D/m ¼ 2pm=N, where m ¼ 0; …; N 1. The case m ¼ 0 corresponds to identical (or in-phase) synchronization considered earlier, while m ¼ 1; 2; …; N 1 correspond to cluster and splay states. The cluster number Nc ¼ lcm ðm; NÞ=m, where lcmð; Þ denotes the least common multiple, determines how many clusters of oscillators exist, with Nc ¼ N corresponding to a splay state.21,59 As it has been shown in Ref. 22, the cluster solution (32) exists if and only if the following conditions hold N X
Gkj cos ½ðj kÞD/m ¼ Gcm ¼ const;
j¼1 N X
Gkj sin ½ðj kÞD/m ¼ Gsm ¼ const;
(33)
j¼1
FIG. 3. Same as Fig. 2 with gamma distributed delay coupling. Parameter values are x ¼ 1, c ¼ 0, h ¼ 0, and l ¼ 1. (a) weak kernel (p ¼ 1), k ¼ 0:1 (black solid), k ¼ 0.25 (blue dashed), k ¼ 0:4 (red dotted), fully synchronized state exists below the corresponding curve; (b) strong kernel (p ¼ 2), k ¼ 0:15 (black solid), k ¼ 0.25 (blue dashed), k ¼ 0:35 (red dotted), fully synchronized state exists below the corresponding curve. (c) and (d) weak kernel (p ¼ 1), k ¼ 0.25, K ¼ 0.5, (c) a ¼ 0:8. (d) a ¼ 1. (e) and (f) strong kernel (p ¼ 2), (e) a ¼ 1:5, (f) a ¼ 3. All eigenvalues of the coupling matrix lie on or inside the black circle.
independently of k. Although these conditions can be satisfied by different kinds of matrices (see Ref. 22 for a detailed discussion), importantly, they are satisfied by all circulant matrices, of which many standard network topologies are special cases. Provided that the above conditions hold, the common radius r0;m and the common frequency Xm are determined by the following equations: 2 ¼ k þ Kl½Gcm ðXm ; hÞ cos h; r0;m
manifold is always unstable for uni-directional coupling (both with and without self-feedback) and is always stable for bi-directional and all-to-all couplings without and with self-feedback. The values of the coupling strength K were taken to be the same as in the earlier work on synchronization of coupled Stuart-Landau oscillators21,22 to illustrate the role played by the width of the delay distribution. Figure 3 illustrates the computation of regions of existence of the fully synchronized solution, and the master stability function for the gamma delay distribution in the case of p ¼ 1 and p ¼ 2. Similarly to the case of the uniform delay distribution, increasing the bifurcation parameter k leads to an increase in the region of parameters, where the fully synchronized solutions exist. For both weak (p ¼ 1) and strong (p ¼ 2) gamma distribution kernels, decreasing the average time delay (which is equivalent to increasing the parameter a in the gamma distribution) leads to a wider region of stability. Unlike the case of the uniform distribution kernel, varying a can lead to (de)stabilization of certain coupling topologies: while for small values of a only the bidirectional and all-to-all couplings can be stable, as a increases the uni-directional coupling is also stabilized. V. STABILITY OF CLUSTER AND SPLAY STATES
So far, we have considered identical or zero-lag synchronization, where all oscillators have exactly the same amplitudes and phases. It is possible to use a similar approach for
(34)
2 Xm ¼ x cr0;m þ Kl½Gsm ðXm ; hÞ sin h;
where Gcm ðXm ; hÞ ¼ Gcm Fc ðXm ; hÞ Gsm Fs ðXm ; hÞ; Gsm ðXm ; hÞ ¼ Gsm Fc ðXm ; hÞ þ Gcm Fs ðXm ; hÞ: Using the ansatz rk ðtÞ ¼ r0;m ½1 þ drk ðtÞ; /k ðtÞ ¼ Xm t þ kD/m þ d/k ðtÞ, the variational equation for linearization near the cluster state (32) can be found as follows: n_ k ¼ ðJ0;m KlGÞnk þ K Gkj
ð1
N X j¼1
gðt0 ÞRj;m ðt0 Þnj ðt t0 Þdt0 ; k ¼ 1; …; N;
(35)
0
where nk ¼ ðdrk ; d/k ÞT , and we have also introduced the matrices 0 J0;m ¼ @
G¼
0
Rj;m ðt Þ ¼
2 2r0;m
0
2 2cr0;m
0
1 A;
Gcm ðXm ; hÞ Gsm ðXm ; hÞ Gsm ðXm ; hÞ cos Um j sin Um j
Gcm ðXm ; hÞ ! sin Um j ; cos Um j
! ;
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 157.211.3.15 On: Mon, 12 Jan 2015 02:35:02
043117-8
€ll Kyrychko, Blyuss, and Scho
Chaos 24, 043117 (2014)
0 and the phase difference Um j ¼ Xm t þ ðj kÞD/m þ h. Using the same approach as before, one can recast this equation in the form ð1 _n ¼ IN ðJ0;m KlGÞn þ K gðt0 Þ½G Rj;m n ðt t0 Þdt0 ; 0
(36) with n ¼ ðn1 ; …; nN ÞT . To make further progress in the analysis of stability of the cluster state, one has to diagonalize the matrix G and study Floquet exponents corresponding to transverse eigenvalues of this matrix. However, unlike the case of identical synchronization m ¼ 0, this is, in general, impossible, unless the matrix G has certain properties which make the matrix Rj;m independent of k. This is the case for all topologies considered in Sec. IV, which are represented by circulant matrices. To illustrate the computation of the master stability function for cluster states, we consider the network of the StuartLandau oscillators with uni-directional ring coupling: Gk;kþ1 ¼ GN;1 ¼ 1 ¼ l, and all other Gkj ¼ 0. In this case, we 0 have Um j ¼ Xm t þ D/m þ h, and the matrix Rj;m reduces to Rj;m ðt0 Þ ¼ R1 ðt0 Þ ¼
cosðh Xm t0 þ D/m Þ sinðh Xm t0 þ D/m Þ sinðh Xm t0 þ D/m Þ
cosðh Xm t0 þ D/m Þ
! :
Since this matrix does not depend on k, one can use the same approach as in the case of a fully synchronized solution to diagonalize the matrix G, which results in the blockdiagonalized variational equation ð1 f_ k ðtÞ ¼ ðJ0;m KlGÞfk ðtÞ þ Kk gðt0 ÞR1 ðt0 Þfðt t0 Þdt0 ;
In the case of a single discrete time delay gðuÞ ¼ dðu sÞ, this equation reduces to a case investigated by Choe et al.21,22 Figure 4 illustrates the computation of the master stability function of the cluster state for uni-directional ring coupling and uniform or gamma delay distributions. For both of these delay distributions, an increase in the number of oscillators appears to increase the range of coupling strengths, for which the cluster state exists, and also proportionally increase the region of stability of this state. In the case of the gamma distributed delay kernel, one can note that whenever it exists, the cluster state is stable for lower values of the coupling strength across the full range of admissible values of parameter a. VI. DISCUSSION
In this paper, we have shown how the framework of the master stability function for the computation of stability of zero-lag and cluster synchronized states in systems of coupled oscillators can be generalized to the case of distributed delay coupling. To illustrate how the master stability function framework can be used for studying the stability of synchronized states, we have considered the example of a network of coupled Stuart-Landau oscillators, which is a normal form of a supercritical Hopf bifurcation. In this example, computation of the master stability function, i.e., the Floquet exponents, reduces to finding the corresponding eigenvalues of a characteristic equation, which can be done semi-analytically in terms of system parameters,
0
where the eigenvalues k of the matrix G are explicitly given by k ¼ e2ikp=N ; k ¼ 0; 1; …; N 1. Since all terms in the above equation are independent of time, the Floquet exponents can be found as the eigenvalues K of the characteristic equation ð1 0 0 Kt0 0 gðt ÞR1 ðt Þe dt ¼ 0; det J0;m KlG KI2 þ Kk 0
(37) where I2 is a 2 2 identity matrix. More explicitly, this transcendental equation has the form 2 þ KðlGcm ðXm ; hÞ k FLc ðXm ; h þ D/m ; KÞÞK K2 þ 2½r0;m 2 K½lGcm ðXm ; hÞ k FLc ðXm ; h þ D/m ; KÞ þ 2r0;m
þ cðlGsm ðXm ; hÞ k FLs ðXm ; h þ D/m ; KÞÞ þ K 2 ½lGcm ðXm ; hÞ k FLc ðXm ; h þ D/m ; KÞ2 þ K 2 ½lGsm ðXm ; hÞ k FLs ðXm ; h þ D/m ; KÞ2 þ K 2 ½lGsm ðXm ; hÞ k FLs ðXm ; h þ D/m ; KÞ2 ¼ 0: (38)
FIG. 4. Master stability function for a splay state with m ¼ 1 for a ring of uni-directionally coupled Stuart-Landau oscillators. Parameter values are x ¼ 1, c ¼ 0, h ¼ 0, and l ¼ 1. (a) and (b) uniform delay distribution, s ¼ 2p; k ¼ 0:1, (a) N ¼ 4, (b) N ¼ 10. (c) and (d) weak kernel (p ¼ 1), k ¼ 0.25, (c) N ¼ 4, (d) N ¼ 10. (e) and (f) strong kernel (p ¼ 2), k ¼ 0.25, (e) N ¼ 4, (f) N ¼ 10. The splay state does not exist in the white region.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 157.211.3.15 On: Mon, 12 Jan 2015 02:35:02
043117-9
€ll Kyrychko, Blyuss, and Scho
topological eigenvalues of the coupling matrix, and the coupling parameters, where the Laplace transform of the delay kernel enters. Calculations for the case of uniform and gamma distributed delay kernels show that increasing the width of the delay distribution (for a uniform kernel) or reducing the mean time delay (for a gamma distributed delay kernel) leads often to a larger region of stability, thus allowing the fully synchronized state to be stable for a larger class of coupling topologies. Although the analysis presented in this paper has focused upon the case of linear coupling between different oscillators in the network, it is straightforward to generalize our results on the stability of zero-lag and cluster synchronization to systems with nonlinear coupling,60 since the master stability equation is still a linear variational equation involving the appropriate delay distribution kernel. Recent studies of amplitude death in systems of delay-coupled oscillators47,48,61 have highlighted the important role played by the coupling phase in determining the regions of possible suppression of oscillations. Due to substantial analytical headway provided by considering Stuart-Landau oscillators, it would be insightful to use the master stability framework developed in this paper to investigate the effects of the coupling phase on existence and stability of synchronized solutions in networks of such oscillators. ACKNOWLEDGMENTS
This work was partially supported by DFG in the framework of SFB 910: Control of self-organizing nonlinear systems: Theoretical methods and concepts of application. Y.N.K. and K.B.B. gratefully acknowledge the hospitality of the Institut f€ ur Theoretische Physik, TU Berlin, where part of this work was completed. E.S. is grateful to Carolin Wille for helpful discussions. The authors thank the anonymous referees for helpful comments and suggestions that have helped to improve the presentation.
1
R. Albert and A.-L. Barabasi, Rev. Mod. Phys. 74, 47 (2002). A.-L. Barabasi and E. Bonabeau, Sci. Am. 288, 60 (2003). 3 M. E. J. Newman, SIAM Rev. 45, 167 (2003). 4 D. J. Watts and S. H. Strogatz, Nature 393, 440 (1998). 5 S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D. U. Hwang, Phys. Rep. 424, 175 (2006). 6 A. Pikovsky, M. Rosenblum, and J. Kurths, Synchronization: A Universal Concept in Nonlinear Sciences (Cambridge University Press, Cambridge, 2001). 7 L. M. Pecora and T. L. Carroll, Phys. Rev. Lett. 80, 2109 (1998). 8 M. G. Earl and S. H. Strogatz, Phys. Rev. E 67, 036204 (2003). 9 O. D’Huys, R. Vicente, T. Erneux, J. Danckaert, and I. Fischer, Chaos 18, 037116 (2008). 10 W. Kinzel, A. Englert, G. Reents, M. Zigzag, and I. Kanter, Phys. Rev. E 79, 056207 (2009). 11 O. D’Huys, R. Vicente, J. Danckaert, and I. Fischer, Chaos 20, 043127 (2010). 12 A. Englert, W. Kinzel, Y. Aviad, M. Butkovski, I. Reidler, M. Zigzag, I. Kanter, and M. Rosenbluh, Phys. Rev. Lett. 104, 114102 (2010). 13 V. Flunkert, S. Yanchuk, T. Dahms, and E. Sch€ oll, Phys. Rev. Lett. 105, 254101 (2010). 14 Complex Time-Delay Systems, edited by F. M. Atay (Springer, New York, 2010). 15 A. Englert, S. Heiligenthal, W. Kinzel, and I. Kanter, Phys. Rev. E 83, 046222 (2011). 2
Chaos 24, 043117 (2014) 16
K. Hicke, O. D’Huys, V. Flunkert, E. Sch€ oll, I. Danckaert, and I. Fischer, Phys. Rev. E 83, 056211 (2011). J. Lehnert, T. Dahms, P. H€ ovel, and E. Sch€ oll, Europhys. Lett. 96, 60013 (2011). 18 V. Flunkert, I. Fischer, and E. Sch€ oll, “Dynamics, control and information in delay-coupled systems,” Philos. Trans. R. Soc. London, Ser. A 371, 20120465 (2013). 19 M. C. Soriano, J. Garcıa-Ojalvo, C. R. Mirasso, and I. Fischer, Rev. Mod. Phys. 85, 421 (2013). 20 A. Zakharova, I. Schneider, Y. N. Kyrychko, K. B. Blyuss, A. Koseska, B. Fiedler, and E. Sch€ oll, Europhys. Lett. 104, 50004 (2013). 21 C.-U. Choe, T. Dahms, P. H€ ovel, and E. Sch€ oll, Phys. Rev. E 81, 025205(R) (2010). 22 C.-U. Choe, T. Dahms, P. H€ ovel, and E. Sch€ oll, in Proceedings of the Eighth AIMS International Conference on Dynamical Systems, Differential Equations and Applications (American Institute of Mathematical Sciences, Springfield, 2011), pp. 292–301. 23 C.-U. Choe, H. Jang, V. Flunkert, T. Dahms, P. H€ ovel, and E. Sch€ oll, Dyn. Syst. 28, 15 (2013). 24 T. Dahms, J. Lehnert, and E. Sch€ oll, Phys. Rev. E 86, 016202 (2012). 25 M. Dhamala, V. K. Jirsa, and M. Ding, Phys. Rev. Lett. 92, 074104 (2004). 26 C. R. S. Williams, T. E. Murphy, R. Roy, F. Sorrentino, T. Dahms, and E. Sch€ oll, Phys. Rev. Lett. 110, 064104 (2013). 27 D. Hunt, G. Korniss, and B. K. Szymanski, Phys. Rev. Lett. 105, 068701 (2010). 28 D. Hunt, B. K. Szymanski, and G. Korniss, Phys. Rev. E 86, 056114 (2012). 29 A. Gjurchinovski and V. Urumov, Europhys. Lett. 84, 40013 (2008). 30 A. Gjurchinovski and V. Urumov, Phys. Rev. E 81, 016209 (2010). 31 A. Gjurchinovski, A. Zakharova, and E. Sch€ oll, Phys. Rev. E 89, 032915 (2014). 32 F. Atay, Phys. Rev. Lett. 91, 094101 (2003). 33 S. A. Campbell and R. Jessop, Math. Model. Nat. Phenom. 4, 1 (2009). 34 C. U. Choe, R.-S. Kim, H. Jang, P. H€ ovel, and E. Sch€ oll, Int. J. Dyn. Control 2, 2 (2014). 35 G. Kiss and B. Krauskopf, Dyn. Syst. 26, 85 (2011). 36 A. Gjurchinovski, T. J€ ungling, V. Urumov, and E. Sch€ oll, Phys. Rev. E 88, 032912 (2013). 37 W. Michiels, V. Van Assche, and S.-I. Niculescu, IEEE Trans. Autom. Control 50, 493 (2005). 38 K. B. Blyuss and Y. N. Kyrychko, Bull. Math. Biol. 72, 490 (2010). 39 T. Faria and S. Trofimchuk, Nonlinearity 23, 2457 (2010). 40 S. A. Gourley and J. W.-H. So, Proc. R. Soc. Edinburgh 133, 527 (2003). 41 R. Sipahi, F. M. Atay, and S.-I. Niculescu, SIAM J. Appl. Math. 68, 738 (2007). 42 C. W. Eurich, A. Thiel, and L. Fahse, Phys. Rev. Lett. 94, 158104 (2005). 43 R. Vicente, L. L. Gollo, C. R. Mirasso, I. Fischer, and P. Gordon, Proc. Natl. Acad. Sci. U.S.A. 105, 17157 (2008). 44 C. Wille, J. Lehnert, and E. Sch€ oll, Phys. Rev. E 90, 032908 (2014). 45 I.-C. Morarescu, W. Michiels, and M. Jungers, in Proceedings of the American Control Conference, Washington (2013). 46 B. Fiedler, V. Flunkert, M. Georgi, P. H€ ovel, and E. Sch€ oll, Phys. Rev. Lett. 98, 114101 (2007). 47 Y. N. Kyrychko, K. B. Blyuss, and E. Sch€ oll, Eur. Phys. J. B 84, 307 (2011). 48 Y. N. Kyrychko, K. B. Blyuss, and E. Sch€ oll, Philos. Trans. R. Soc. London, Ser. A 371, 20120466 (2013). 49 K. Pyragas and T. Pyragiene, Phys. Rev. E 78, 046217 (2008). 50 M. Adimy, F. Crauste, and S. Ruan, Nonlinear Anal.: Real World Appl. 6, 651 (2005). 51 M. Barrio, K. Burrage, A. Leier, and T. Tian, PLoS Comput. Biol. 2, e117 (2006). 52 S. P. Blythe, R. M. Nisbet, W. S. C. Gurney, and N. MacDonald, J. Math. Anal. Appl. 109, 388 (1985). 53 K. Cooke and Z. Grossman, J. Math. Anal. Appl. 86, 592 (1982). 54 J. M. Cushing, in Mathematics of Biology, edited by M. Iannelli (Springer, New York, 2011). 55 G. Stepan, in Dynamics and Chaos in Manufacturing Process, edited by F. C. Moon (Wiley, New York, 1998). 56 J. E. Mittler, B. Sulzer, A. U. Neumann, and A. S. Perelson, Math. Biosci. 152, 143 (1998). 57 R. Sipahi, S.-I. Niculescu, and F. M. Atay, in Proceedings of the American Control Conference, New York (2007). 17
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 157.211.3.15 On: Mon, 12 Jan 2015 02:35:02
043117-10 58
€ ll Kyrychko, Blyuss, and Scho
O. Roesch, H. Roth, and S.-I. Niculescu, in Proceedings of the IEEE International Conference on Mechatronics and Automation, Niagara Falls (2005). 59 R. Zillmer, R. Livi, A. Politi, and A. Torcini, Phys. Rev. E 76, 046102 (2007).
Chaos 24, 043117 (2014) 60
61
O. V. Popovych, C. Hauptmann, and P. A. Tass, Phys. Rev. Lett. 94, 164102 (2005). W. Zou, J. Lu, Y. Tang, C. Zhang, and J. Kurths, Phys. Rev. E 84, 066208 (2011).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 157.211.3.15 On: Mon, 12 Jan 2015 02:35:02
