PDF - Chaos at Maryland - University of Maryland
PHYSICAL REVIEW E 77, 036107 共2008兲
Synchronization in networks of networks: The onset of coherent collective behavior in systems of interacting populations of heterogeneous oscillators Ernest Barreto,1,* Brian Hunt,2 Edward Ott,3 and Paul So1
1
Department of Physics & Astronomy, The Center for Neural Dynamics and The Krasnow Institute for Advanced Study, George Mason University, Fairfax, Virginia 22030, USA 2 Department of Mathematics and Institute for Physical Science and Technology, University of Maryland, College Park, Maryland 20742, USA 3 Institute for Research in Electronics and Applied Physics, Department of Physics, and Department of Electrical and Computer Engineering, University of Maryland, College Park, Maryland 20742, USA 共Received 7 August 2007; revised manuscript received 10 December 2007; published 6 March 2008兲 The onset of synchronization in networks of networks is investigated. Specifically, we consider networks of interacting phase oscillators in which the set of oscillators is composed of several distinct populations. The oscillators in a given population are heterogeneous in that their natural frequencies are drawn from a given distribution, and each population has its own such distribution. The coupling among the oscillators is global, however, we permit the coupling strengths between the members of different populations to be separately specified. We determine the critical condition for the onset of coherent collective behavior, and develop the illustrative case in which the oscillator frequencies are drawn from a set of 共possibly different兲 Cauchy-Lorentz distributions. One motivation is drawn from neurobiology, in which the collective dynamics of several interacting populations of oscillators 共such as excitatory and inhibitory neurons and glia兲 are of interest. DOI: 10.1103/PhysRevE.77.036107
PACS number共s兲: 89.75.Fb, 05.45.Xt, 05.45.⫺a, 87.19.L⫺
In recent years, there has been considerable interest in networks of interacting systems. Researchers have found that an appropriate description of such systems involves an understanding of both the dynamics of the individual oscillators and the connection topology of the network. Investigators studying the latter have found that many complex networks have a modular structure involving motifs 关1兴, communities 关2,3兴, layers 关4兴, or clusters 关5兴. For example, recent work has shown that as many kinds of networks 共including isotropic homogeneous networks and a class of scale-free networks兲 transition to full synchronization, they pass through epochs in which well-defined synchronized communities appear and interact with one another 关3兴. Knowledge of this structure, and the dynamical behavior it supports, informs our understanding of biological 关6兴, social 关7兴, and technological networks 关8兴. Here we consider the onset of coherent collective behavior in similarly structured systems for which the dynamics of the individual oscillators is simple. In seminal work, Kuramoto analyzed a mathematical model that illuminates the mechanisms by which synchronization arises in a large set of globally coupled phase oscillators 关9兴. An important feature of Kuramoto’s model is that the oscillators are heterogeneous in their frequencies. And, although these mathematical results assume global coupling, they have been fruitfully applied to further our understanding of systems of fireflies, arrays of Josephson junctions, electrochemical oscillators, and many other cases 关10兴. In this work, we study systems of several interacting Kuramoto systems, i.e., networks of interacting populations of phase oscillators. Our motivation is drawn not only from the applications listed above 共e.g., an amusing application
*[email protected] 1539-3755/2008/77共3兲/036107共7兲
might be interacting populations of fireflies, where each population inhabits a separate tree兲, but also from other biological systems. Rhythms arising from coupled cell populations are seen in many of the body’s organs 共including the heart, the pancreas, and the kidneys, to name but a few兲, all of which interact. For example, the circadian rhythm influences many of these systems. We draw additional motivation from consideration of how neuronal tissue is organized. Although we do not consider neuronal systems specifically in this paper, we note that heterogeneous ensembles of neurons often exhibit a “network-of-networks” topology. At the cellular level, populations of particular kinds of neurons 共e.g., excitatory neurons兲 interact not only among themselves, but also with populations of other distinct neuronal types 共e.g., inhibitory neurons兲. At a higher level of organization, various anatomical regions of the brain interact with one another as well 关6兴. Although our network is simple, it is novel in that we include heterogeneity at several levels. Each population consists of a collection of phase oscillators whose natural frequencies are drawn at random from a given distribution. To allow for heterogeneity at the population level, we let each population have its own such frequency distribution. In addition, our system is heterogeneous at the coupling level: we consider global coupling such that the coupling strengths between the members of different populations can be separately specified. The assumption of global 共but population-weighted兲 coupling permits an analytic determination of the critical condition for the onset of coherent collective behavior, as we will show. While this assumption may not strictly apply in some applications, our results provide insight into the behavior of networks of similar topology even if the connectivity is less than global. We begin by specifying our model and deriving our main
036107-1
©2008 The American Physical Society
PHYSICAL REVIEW E 77, 036107 共2008兲
BARRETO et al.
results. We then discuss several illustrative examples. Consider first a two-species Kuramoto model. We label the separate species and and assume that there are N and N such oscillators in each population, respectively. Thus the system equations are
G 共 兲 =
+
and the order parameters are given by r e i =
k 兺 sin共 j − i − ␣兲, N j=1
冉
N
⬘
共3兲
Here, is an overall coupling strength, the ␣’s provide additional heterogeneity in the coupling functions, and the matrix
冊
共1兲
defines the connectivity among the populations 关11兴. For arbitrarily many different species, let range over the various population symbols , , . . . with the understanding that depending on the context, may represent either a label 共when subscripted兲 or a variable. Thus we have
冋
册
k⬘ ⬘ di = i + 兺 兺 sin共⬘j − i − ␣⬘兲 . N⬘ j=1 dt ⬘
The incoherent state has distributed uniformly over 关0 , 2兴, so that f = 1 / 2 and r = 0. These satisfy Eq. 共3兲 trivially. We test the stability of this solution by perturbing it. Note that a perturbation to f leads to a perturbation of r, and we expect these perturbations to either grow or shrink exponentially in time, depending on the overall coupling strength . The marginally stable case occurs at a critical value ⴱ at which coherent collective behavior emerges. Thus we write f = 1 / 2 + 共␦ f 兲est and r = 共␦r兲est, where 共␦ f 兲 and 共␦r兲 are small. Inserting the first of these into the continuity equation 关Eq. 共3兲兴, and keeping only first-order terms, we obtain s共␦ f 兲 +
The i are the constant natural frequencies of the oscillators when uncoupled, and are distributed according to a set of time-independent distribution functions G共兲 关12兴. We define the usual Kuramoto order parameter for each population, i.e., r e i =
冊
f 兵关 + 兺 k⬘r⬘ sin共⬘ − − ␣⬘兲兴f 其 = 0. + t
N
N
共2兲
and if we write F共 , , t兲 = f 共 , , t兲G共兲, we have
k + 兺 sin共 j − i − ␣兲. N j=1
k k k k
F e id d .
F ជ d + ⵜ · F ˆ = 0, dt t
k di = i + 兺 sin共 j − i − ␣兲 N j=1 dt
冉
冕冕
In this context, the distribution functions satisfy the equations of continuity, i.e.,
N
K=
F共, ,t兲d ,
0
N
k di = i + 兺 sin共 j − i − ␣兲 N j=1 dt
冕
2
⬘
⫻cos共⬘ − − ␣⬘兲. The solution to this equation is 共 ␦ f 兲 =
N
1 兺 e i j . N j=1
1 共 ␦ f 兲 = 兺 k⬘共␦r⬘兲 2
冋
册
1 ei共⬘−−␣⬘兲 e−i共⬘−−␣⬘兲 k 共 ␦ r 兲 兺 ⬘ ⬘ s − i + s + i . 4 ⬘
共4兲
Here, r describes the degree of synchronization in each population, and ranges from 0 to 1. Using this, the above equations can be expressed as
Consistency demands that the perturbations 共␦ f 兲 and 共␦r兲 be related to each other via the order parameter equation, Eq. 共2兲. This yields our main result, as follows. Equation 共2兲 becomes
di = i + 兺 k⬘r⬘ sin共⬘ − i − ␣⬘兲. dt
共␦r兲estei =
⬘
冕
⬁
−⬁
Assuming that the N are very large, we solve for the onset of coherent collective behavior by using a distribution function approach. Thus instead of discrete indices, we imagine continua of oscillators described by distribution functions F共 , , t兲 such that F共 , , t兲dd is the fraction of oscillators whose phases and natural frequencies lie in the infinitesimal volume element dd at time t. Note that in the N → ⬁ limit,
G 共 兲
冕冉 2
0
冊
1 + 共␦ f 兲est eidd . 2
The integral involving 1 / 2 is zero. Inserting the solution for 共␦ f 兲 from Eq. 共4兲, one obtains 关13兴 共␦r兲ei =
冋冕 1 2
⬁
−⬁
G 共 兲 d s − i
册兺 ⬘
k⬘共␦r⬘兲ei共⬘−␣⬘兲 .
Define the bracketed expression to be g共s兲, and define z = 共␦r兲ei. Then, we have
036107-2
PHYSICAL REVIEW E 77, 036107 共2008兲
SYNCHRONIZATION IN NETWORKS OF NETWORKS: THE ...
z = g共s兲 兺 k⬘e−i␣⬘z⬘ . ⬘
Now define the complex quantity ¯k⬘ = k⬘e−i␣⬘. Using the Kronecker delta ␦⬘, the above equation can be written
兺 ⬘
冉
¯k − ⬘
␦⬘
g共s兲
冊
tions; we set ⌬ = ⌬ = ⌬ and ⍀ = ⍀ = ⍀. 共A more generic example will follow.兲 Denoting D = det共K兲 and T = tr共K兲, we separate the real and imaginary parts of Eq. 共7兲 to obtain D2 − 2⌬T + 4⌬2 − 4共v + ⍀兲2 = 0, 共v + ⍀兲共4⌬ − T兲 = 0.
z⬘ = 0.
One solution to these equations is
In matrix notation for the case of two populations labeled and , this is
冢
¯k −
1 g共s兲
¯k
¯k ¯k − 1 g共s兲
冣
冉冊
z = 0. z
共5兲
This equation has a nontrivial solution if the determinant of the matrix is zero. This condition determines the growth rate s in terms of , ¯k⬘, and the parameters of the distributions G共兲 关via g共s兲兴. To illustrate the resulting behavior, we note that g共s兲 can be evaluated in closed form for a Cauchy-Lorentz distribution G 共 兲 =
1 ⌬ , 共 − ⍀兲2 + ⌬2
共6兲
where ⍀, the mean frequency of population , and the half width at half maximum ⌬ are both real, and ⌬ is positive. One obtains
ⴱ = ⌬
vⴱ = − ⍀,
Using this, the determinant condition for the two-population case reduces to ¯ − 2共s + ⌬ + i⍀ 兲兴关k ¯ − 2共s + ⌬ + i⍀ 兲兴 − ¯k ¯k = 0. 关k For simplicity, we set the phase angles ␣⬘ to zero for the remainder of this paper 关14兴. In this case, the matrix elements ¯k ¯ ⬘ are purely real, so that k⬘ = k⬘. The determinant condition then becomes 关k − 2共s + ⌬ + i⍀兲兴关k − 2共s + ⌬ + i⍀兲兴 共7兲
Notice that if = 0, then s = −⌬ − i⍀, indicating that the incoherent state is stable for zero coupling 共since −⌬ is negative兲. We imagine increasing 共or decreasing兲 until s crosses the imaginary axis at a critical value ⴱ. At this point, the incoherent state loses stability and coherent collective behavior emerges in the ensemble. Thus the critical value ⴱ can be determined from the determinant condition by setting s = iv, where v is real 共so that the perturbations are marginally stable兲, and equating the real and imaginary parts of the left side of Eq. 共7兲 to zero. This results in two equations which can be solved simultaneously for the two 共real兲 unknowns and v. For our first example, we choose two identical popula-
冉
冊
T ⫾ 冑T2 − 4D , D
which is valid for T2 ⱖ 4D, since is assumed to be real. Notice that the appropriate solution as D → 0 共using the negative sign for T ⬎ 0 and the positive sign for T ⬍ 0兲 is vⴱ = − ⍀,
ⴱ =
2⌬ , T
as can be verified by setting D = 0 in Eq. 共8兲 directly. Another solution is vⴱ = − ⍀ ⫾
⌬ 冑4D − T2, T
ⴱ =
4⌬ . T
Finally, setting T = 0 in Eq. 共8兲 yields vⴱ = −⍀ and ⴱ ⌬ = ⫾ 冑2−D for D ⬍ 0, and no solution for D ⱖ 0. These results are summarized in Table I. Thus the critical values ⴱ are determined by T, D, and ⌬. To illustrate this result, we begin by discussing a particular example. Consider the matrix K=
1 = 2共s + ⌬ + i⍀兲. g 共 兲
− 2kk = 0.
共8兲
冉 冊 1 −1 1
0
,
which has trace T = 1 and determinant D = 1. This corresponds to case 2 in Table I, from which we find that ⴱ = 4⌬ / T. Figure 1 shows the results of a numerical simulation of two populations of 10 000 oscillators each, using ⌬ = 1. The order parameters r vs are shown, and we can see that the oscillator populations are incoherent for values below the predicted critical value ⴱ = 4, and that they grow increasingly synchronized as is increased beyond ⴱ. Next, we examined eight different connectivity matrices K that were chosen to sample the various regions in T − D space. Table II shows these matrices and the corresponding value共s兲 of ⴱ. These predictions were tested by numerically calculating the order parameters r for a range of coupling TABLE I. Solutions to Eq. 共8兲 for two identical populations. D = det共K兲 and T = tr共K兲, where K is the connectivity matrix 关Eq. 共1兲兴, and ⌬ is the width parameter in Eq. 共6兲. Case
Condition
vⴱ
1
T2 ⬎ 4D
−⍀
2
T ⱕ 4D
3
D=0
036107-3
2
−⍀ ⫾ T 冑4D − T −⍀ ⌬
4
T = 0, D ⬍ 0
−⍀
5
T = 0, D ⱖ 0
no solution
ⴱ ⌬共 2
T⫾冑T2−4D
D
兲
4⌬
T
2⌬
T
2⌬
⫾ 冑−D no solution
PHYSICAL REVIEW E 77, 036107 共2008兲
BARRETO et al.
TABLE II. Connectivity matrices K chosen to sample T-D space. Case A
B
C
D FIG. 1. Numerical calculation of the order parameters 共쎲 , △兲 vs for case A in Table II. The vertical line corresponds to the predicted value ⴱ = 4. The data point nearest ⴱ is at = 4.15.
values 关15兴. Results are shown in Fig. 2. The various cases are located by letter in the T − D plane according to the trace and determinant of the matrices, and the corresponding inset shows the numerically calculated order parameters plotted vs , with the predicted critical coupling indicated by a vertical line at = ⴱ. For example, the inset to case A 关with 共T , D兲 = 共1 , 1兲兴 corresponds to Fig. 1, which was discussed above. It can be seen that for all cases, the onset of synchronization occurs as predicted. Note that more than one prediction for ⴱ may be specified by our analysis 共see Table I兲. The solutions closest to = 0 are the relevant ones, because we expect the incoherent state to lose stability once the first ⴱ solution is encountered. There are two possible cases depending on the sign of D. First, if the two solutions have the same sign, then there is only one critical value ⴱ 共which may be positive or negative depending on the sign of the trace兲 for the onset of synchronization. This occurs for D ⬎ 0 and T ⫽ 0, as can be seen in
E
F
G
H
T
D
ⴱ
1 −1 1 0
1
1
4
−2 −3 1 1
−1
1
−4
3 1 −3.5 −1
2
0.5
2共2 − 冑2兲 = 1.172
−3 1 −3.5 1
−2
0.5
2共−2 + 冑2兲 = −1.172
−1 −1 1 2
1
−1
−共1 ⫾ 冑5兲 = −3.236, 1.236
1 1 −1 −2
−1
−1
1 ⫾ 冑5 = −1.236, 3.236
2 1 −3 −2
0
−1
⫾2
1 −1 2 −1
0
1
none
Matrix
共 兲 共 兲
共 兲 共 兲 共 兲 共 兲 共 兲 共 兲
Fig. 2. 共Interestingly, for D ⬎ 0 and T = 0, synchronization does not occur for any .兲 The other case occurs for D ⬍ 0, for which the two ⴱ solutions have opposite signs. In this case, there are two critical values ⴱ for the onset of synchronization—one on either side of = 0. This can also be observed in Fig. 2. In the more general case in which the various populations have different natural frequency distributions, it is not typi-
FIG. 2. Numerical simulations using the matrices listed in Table II for identical populations. The letters indicate the placement of each case in the T-D plane, and the corresponding insets show numerical calculations of the order parameters 共쎲 , △兲 vs for that case 共in all cases, = 0 is in the center of the horizontal axis兲. Vertical lines in the insets indicate the predicted value共s兲 ⴱ listed in Table II for the onset of coherent collective behavior. In all cases, we used ⌬ = 1. Note that for D ⬎ 0, there is one value of ⴱ, whose sign corresponds to the sign of the trace T. If D ⱖ 0 and T = 0, then synchronization is not possible for any coupling strength. For D ⬍ 0, there are two values of ⴱ: one positive, and one negative.
036107-4
PHYSICAL REVIEW E 77, 036107 共2008兲
SYNCHRONIZATION IN NETWORKS OF NETWORKS: THE ... 6
6
4
Im(η)
Im(η)
4 2 0 −2
−4 −4
−3
v
−2
−1
−5
−4
−6
−5
−4
v
−3
−2
−1
−3
−2
−1
6
0
4
−1
Re(η)
Re(η)
−6
(a)
1
−2
2 0
−3 −4
0 −2
(a)
(b)
2
−4
−3
v
−2
−2
−1
(b)
FIG. 3. Case E, different populations. The upper panel shows Im共共1,2兲兲, with the vertical lines identifying roots at v1 = −4.024 and v2 = −1.722. The lower panel shows Re共共1,2兲兲; values at the roots found above are indicated by horizontal lines, yielding 共1兲 ⴱ = 0.515 and 共2兲 ⴱ = −2.809. Thus we expect synchronization to occur at these values as is either increased or decreased away from zero.
cally possible to describe the onset of synchronization in terms of the determinant and trace of the coupling matrix K = k⬘ alone. We now consider this situation, but retain the Cauchy-Lorentz form of the natural frequency distributions for convenience. We manipulate Eq. 共7兲 as follows. Let s = iv 共i.e., purely imaginary, to consider the marginally stable case兲 and define a = ⌬ + i共v + ⍀兲 and b = ⌬ + i共v + ⍀兲. We obtain
v
FIG. 5. Case A, different populations. Panels as in Fig. 3. Note that the standard branch cut leads to discontinuities and the occurrence of two roots for Im共共1兲兲共v兲 共solid lines, upper panel兲, and none for Im共共2兲兲共v兲 共dotted lines, upper panel兲. From the lower panel we find ⴱ = 2.189, 4.501, taking care to obtain these from Re共共1兲兲 共solid line, lower panel兲. We expect to find synchronization onset at the smaller of these values.
There are two unknowns in this equation: and v. Equation 共9兲 is a quadratic equation in with complex coefficients,
and we can easily obtain two complex solutions 共1,2兲 as functions of v. Since the critical values 共1,2兲 must be real, ⴱ we solve for the roots of Im共共1,2兲兲, and evaluate Re共共1,2兲兲 at . these roots. This yields the possible critical values 共1,2兲 ⴱ Typically, these steps must be performed symbolically and/or numerically; we used MATLAB® 关16兴. As before, the values of 共1,2兲 that are closest to zero 共on either side兲 are the relⴱ evant values. To illustrate, we choose two populations with CauchyLorentz natural frequency distributions 关Eq. 共6兲兴 with parameters ⌬ = 1, ⍀ = 2, ⌬ = 0.5, and ⍀ = 4. We consider the
FIG. 4. Case E, different populations. Calculations of the order parameters vs confirm that the onset of synchronization occurs at ⴱ = −2.809 and 0.515 共vertical lines兲, as predicted in Fig. 3.
FIG. 6. Case A, different populations: The onset of synchronization occurs at ⴱ = 2.189 共vertical line兲, as predicted in Fig. 5. No synchronization is observed for smaller values of 共not shown兲.
D2 − 2共bk + ak兲 + 4ab = 0.
共9兲
036107-5
PHYSICAL REVIEW E 77, 036107 共2008兲
6
3
4
2
Im(η)
Im(η)
BARRETO et al.
2 0
−6 (a)
1 0 −1
−5
−4
−3
v
−2
−1
0
−2
1 (a)
−4
−3
−4
−3
−2
−1
0
−2
−1
0
v
6
3
4
2
2
Re(η)
Re(η)
8
0 −2 −4
(b)
−6
FIG. 7. Case = −1.429, 5.000.
−5
H,
−4
−3
v
different
−2
−1
0
populations.
0 −1
1
We
1
find
ⴱ
same K matrices as before, i.e., those listed in the second column of Table II. Case E is straightforward; the analysis is illustrated in Fig. 3 and the numerical verification is in Fig. 4. Note that since Eq. 共9兲 has complex coefficients, obtaining typically involves taking the square root of a complex number. Therefore one must be mindful of branch cuts when obtaining symbolic and/or numerical solutions. This is important in the analysis for case A, shown in Figs. 5 and 6. Finally, case H, which exhibits no synchronization for identical populations for any value of , does show synchronization in the present case. The analysis is shown in Figs. 7 and 8. We close by giving an example with three different populations. We choose the same Cauchy-Lorentz distributions as above, and add a third with ⌬ = 1 / 3 and ⍀ = 1. We use the following K matrix:
FIG. 8. Case H, different populations. Synchronization occurs at
ⴱ = −1.429 and 5.000 共vertical lines兲, as predicted in Fig. 7.
(b)
v
FIG. 9. 共Color online兲 Three populations. The imaginary and real parts of 共1,2,3兲 are plotted to illustrate the discontinuities due to branch cuts. The analysis proceeds as in the previous cases. Because of branch cuts, two roots occur on 共1兲, one on 共2兲, and none on 共3兲. Evaluating the real parts appropriately, we find ⴱ = −0.891, −0.564, 2.303.
冢
−1
1
1
1
−1
1
1
1
−1
冣
.
The procedure for deriving ⴱ proceeds as above, except that Eq. 共7兲 is replaced by a third-degree polynomial in . Figure 9 shows the imaginary and real parts of the three solutions. 共Note that the branch cuts are more complicated.兲 The predicted onset of synchronization was verified, as shown in Fig. 10.
FIG. 10. Three populations. Synchronization occurs at ⴱ = −0.564 and 2.303, as predicted in Fig. 9.
036107-6
PHYSICAL REVIEW E 77, 036107 共2008兲
SYNCHRONIZATION IN NETWORKS OF NETWORKS: THE ...
ACKNOWLEDGMENTS
In conclusion, we have described how to determine the onset of coherent collective behavior in systems of interacting Kuramoto systems, i.e., systems of interacting populations of phase oscillators with both node and coupling heterogeneity.
E.B. was supported by NIH Grant No. R01-MH79502; E.O. was supported by ONR 共Physics兲 and by NSF Grant No. PHY045624.
关1兴 R. Milo et al., Science 298, 824 共2002兲. 关2兴 M. E. J. Newman and M. Girvan, Phys. Rev. E 69, 026113 共2004兲. 关3兴 A. Arenas, A. Díaz-Guilera, and C. J. Pérez-Vicente, Phys. Rev. Lett. 96, 114102 共2006兲; Physica D 224, 27 共2006兲. 关4兴 M. Kurant and P. Thiran, Phys. Rev. Lett. 96, 138701 共2006兲. 关5兴 X. Wang, L. Huang, Y.-C. Lai, and C. H. Lai, Phys. Rev. E 76, 056113 共2007兲. 关6兴 C. Zhou, L. Zemanová, G. Zamora, C. C. Hilgetag, and J. Kurths, Phys. Rev. Lett. 97, 238103 共2006兲. 关7兴 M. Girvan and M. E. J. Newman, Proc. Natl. Acad. Sci. U.S.A. 99, 7821 共2002兲. 关8兴 R. Milo et al., Science 303, 1538 共2004兲. 关9兴 Y. Kuramoto, in International Symposium on Mathematical Problems in Theoretical Physics, edited by H. Araki, Lecture Notes in Physics Vol. 39 共Springer, New York, 1975兲, 420; Chemical Oscillators, Waves and Turbulence 共Springer, Berlin, 1984兲. 关10兴 A. T. Winfree, The Geometry of Biological Time 共Springer, New York, 1980兲; S. H. Strogatz, Physica D 143, 1 共2000兲; J. A. Acebrón et al.; Rev. Mod. Phys. 77, 137 共2005兲; K. Wiesenfeld and J. W. Swift, Phys. Rev. E 51, 1020 共1995兲; I. Z. Kiss, Y. Zhai, and J. L. Hudson, Science 296, 1676 共2002兲. 关11兴 A similar system of two asymmetrically interacting populations with a particular form of coupling matrix K was consid-
ered in E. Montbrio, J. Kurths, and B. Blasius, Phys. Rev. E 70, 056125 共2004兲. The formulation in the current paper is more general in two important aspects: K is arbitrary, and we allow for any number of interacting populations. Our system is similar to that studied in J. G. Restrepo, E. Ott, and B. R. Hunt, Chaos 16, 015107 共2006兲. However, our formulation permits different natural frequency distributions for each population. The bracketed expression is valid for Re共s兲 ⬎ 0 and may be analytically continued into the region where Re共s兲 ⱕ 0. See E. Ott, P. So, E. Barreto, and T. Antonsen, Physica D 173, 29 共2002兲; E. Ott, Chaos in Dynamical Systems, 2nd ed. 共Cambridge University Press, Cambridge, England, 2002兲, p. 240. Many interesting states require ␣⬘ ⫽ 0, such as the chimera state observed in D. M. Abrams and S.H. Strogatz, Int. J. Bifurcation Chaos Appl. Sci. Eng. 16, 21 共2006兲. Simulations used fourth-order Runge-Kutta with a time step of 0.01 s, N = 10 000 or 50 000, and parameters as noted. The system was initialized in the incoherent state and an initial transient was discarded. The order parameters were then averaged over the subsequent 10 s. Because the standard deviations over this interval were small, no error bars were plotted. MATLAB®, The MathWorks, Inc., Natick, MA, http:// www.mathworks.com/
关12兴
关13兴
关14兴
关15兴
关16兴
036107-7
Synchronization in networks of networks: The onset of coherent collective behavior in systems of interacting populations of heterogeneous oscillators Ernest Barreto,1,* Brian Hunt,2 Edward Ott,3 and Paul So1
1
Department of Physics & Astronomy, The Center for Neural Dynamics and The Krasnow Institute for Advanced Study, George Mason University, Fairfax, Virginia 22030, USA 2 Department of Mathematics and Institute for Physical Science and Technology, University of Maryland, College Park, Maryland 20742, USA 3 Institute for Research in Electronics and Applied Physics, Department of Physics, and Department of Electrical and Computer Engineering, University of Maryland, College Park, Maryland 20742, USA 共Received 7 August 2007; revised manuscript received 10 December 2007; published 6 March 2008兲 The onset of synchronization in networks of networks is investigated. Specifically, we consider networks of interacting phase oscillators in which the set of oscillators is composed of several distinct populations. The oscillators in a given population are heterogeneous in that their natural frequencies are drawn from a given distribution, and each population has its own such distribution. The coupling among the oscillators is global, however, we permit the coupling strengths between the members of different populations to be separately specified. We determine the critical condition for the onset of coherent collective behavior, and develop the illustrative case in which the oscillator frequencies are drawn from a set of 共possibly different兲 Cauchy-Lorentz distributions. One motivation is drawn from neurobiology, in which the collective dynamics of several interacting populations of oscillators 共such as excitatory and inhibitory neurons and glia兲 are of interest. DOI: 10.1103/PhysRevE.77.036107
PACS number共s兲: 89.75.Fb, 05.45.Xt, 05.45.⫺a, 87.19.L⫺
In recent years, there has been considerable interest in networks of interacting systems. Researchers have found that an appropriate description of such systems involves an understanding of both the dynamics of the individual oscillators and the connection topology of the network. Investigators studying the latter have found that many complex networks have a modular structure involving motifs 关1兴, communities 关2,3兴, layers 关4兴, or clusters 关5兴. For example, recent work has shown that as many kinds of networks 共including isotropic homogeneous networks and a class of scale-free networks兲 transition to full synchronization, they pass through epochs in which well-defined synchronized communities appear and interact with one another 关3兴. Knowledge of this structure, and the dynamical behavior it supports, informs our understanding of biological 关6兴, social 关7兴, and technological networks 关8兴. Here we consider the onset of coherent collective behavior in similarly structured systems for which the dynamics of the individual oscillators is simple. In seminal work, Kuramoto analyzed a mathematical model that illuminates the mechanisms by which synchronization arises in a large set of globally coupled phase oscillators 关9兴. An important feature of Kuramoto’s model is that the oscillators are heterogeneous in their frequencies. And, although these mathematical results assume global coupling, they have been fruitfully applied to further our understanding of systems of fireflies, arrays of Josephson junctions, electrochemical oscillators, and many other cases 关10兴. In this work, we study systems of several interacting Kuramoto systems, i.e., networks of interacting populations of phase oscillators. Our motivation is drawn not only from the applications listed above 共e.g., an amusing application
*[email protected] 1539-3755/2008/77共3兲/036107共7兲
might be interacting populations of fireflies, where each population inhabits a separate tree兲, but also from other biological systems. Rhythms arising from coupled cell populations are seen in many of the body’s organs 共including the heart, the pancreas, and the kidneys, to name but a few兲, all of which interact. For example, the circadian rhythm influences many of these systems. We draw additional motivation from consideration of how neuronal tissue is organized. Although we do not consider neuronal systems specifically in this paper, we note that heterogeneous ensembles of neurons often exhibit a “network-of-networks” topology. At the cellular level, populations of particular kinds of neurons 共e.g., excitatory neurons兲 interact not only among themselves, but also with populations of other distinct neuronal types 共e.g., inhibitory neurons兲. At a higher level of organization, various anatomical regions of the brain interact with one another as well 关6兴. Although our network is simple, it is novel in that we include heterogeneity at several levels. Each population consists of a collection of phase oscillators whose natural frequencies are drawn at random from a given distribution. To allow for heterogeneity at the population level, we let each population have its own such frequency distribution. In addition, our system is heterogeneous at the coupling level: we consider global coupling such that the coupling strengths between the members of different populations can be separately specified. The assumption of global 共but population-weighted兲 coupling permits an analytic determination of the critical condition for the onset of coherent collective behavior, as we will show. While this assumption may not strictly apply in some applications, our results provide insight into the behavior of networks of similar topology even if the connectivity is less than global. We begin by specifying our model and deriving our main
036107-1
©2008 The American Physical Society
PHYSICAL REVIEW E 77, 036107 共2008兲
BARRETO et al.
results. We then discuss several illustrative examples. Consider first a two-species Kuramoto model. We label the separate species and and assume that there are N and N such oscillators in each population, respectively. Thus the system equations are
G 共 兲 =
+
and the order parameters are given by r e i =
k 兺 sin共 j − i − ␣兲, N j=1
冉
N
⬘
共3兲
Here, is an overall coupling strength, the ␣’s provide additional heterogeneity in the coupling functions, and the matrix
冊
共1兲
defines the connectivity among the populations 关11兴. For arbitrarily many different species, let range over the various population symbols , , . . . with the understanding that depending on the context, may represent either a label 共when subscripted兲 or a variable. Thus we have
冋
册
k⬘ ⬘ di = i + 兺 兺 sin共⬘j − i − ␣⬘兲 . N⬘ j=1 dt ⬘
The incoherent state has distributed uniformly over 关0 , 2兴, so that f = 1 / 2 and r = 0. These satisfy Eq. 共3兲 trivially. We test the stability of this solution by perturbing it. Note that a perturbation to f leads to a perturbation of r, and we expect these perturbations to either grow or shrink exponentially in time, depending on the overall coupling strength . The marginally stable case occurs at a critical value ⴱ at which coherent collective behavior emerges. Thus we write f = 1 / 2 + 共␦ f 兲est and r = 共␦r兲est, where 共␦ f 兲 and 共␦r兲 are small. Inserting the first of these into the continuity equation 关Eq. 共3兲兴, and keeping only first-order terms, we obtain s共␦ f 兲 +
The i are the constant natural frequencies of the oscillators when uncoupled, and are distributed according to a set of time-independent distribution functions G共兲 关12兴. We define the usual Kuramoto order parameter for each population, i.e., r e i =
冊
f 兵关 + 兺 k⬘r⬘ sin共⬘ − − ␣⬘兲兴f 其 = 0. + t
N
N
共2兲
and if we write F共 , , t兲 = f 共 , , t兲G共兲, we have
k + 兺 sin共 j − i − ␣兲. N j=1
k k k k
F e id d .
F ជ d + ⵜ · F ˆ = 0, dt t
k di = i + 兺 sin共 j − i − ␣兲 N j=1 dt
冉
冕冕
In this context, the distribution functions satisfy the equations of continuity, i.e.,
N
K=
F共, ,t兲d ,
0
N
k di = i + 兺 sin共 j − i − ␣兲 N j=1 dt
冕
2
⬘
⫻cos共⬘ − − ␣⬘兲. The solution to this equation is 共 ␦ f 兲 =
N
1 兺 e i j . N j=1
1 共 ␦ f 兲 = 兺 k⬘共␦r⬘兲 2
冋
册
1 ei共⬘−−␣⬘兲 e−i共⬘−−␣⬘兲 k 共 ␦ r 兲 兺 ⬘ ⬘ s − i + s + i . 4 ⬘
共4兲
Here, r describes the degree of synchronization in each population, and ranges from 0 to 1. Using this, the above equations can be expressed as
Consistency demands that the perturbations 共␦ f 兲 and 共␦r兲 be related to each other via the order parameter equation, Eq. 共2兲. This yields our main result, as follows. Equation 共2兲 becomes
di = i + 兺 k⬘r⬘ sin共⬘ − i − ␣⬘兲. dt
共␦r兲estei =
⬘
冕
⬁
−⬁
Assuming that the N are very large, we solve for the onset of coherent collective behavior by using a distribution function approach. Thus instead of discrete indices, we imagine continua of oscillators described by distribution functions F共 , , t兲 such that F共 , , t兲dd is the fraction of oscillators whose phases and natural frequencies lie in the infinitesimal volume element dd at time t. Note that in the N → ⬁ limit,
G 共 兲
冕冉 2
0
冊
1 + 共␦ f 兲est eidd . 2
The integral involving 1 / 2 is zero. Inserting the solution for 共␦ f 兲 from Eq. 共4兲, one obtains 关13兴 共␦r兲ei =
冋冕 1 2
⬁
−⬁
G 共 兲 d s − i
册兺 ⬘
k⬘共␦r⬘兲ei共⬘−␣⬘兲 .
Define the bracketed expression to be g共s兲, and define z = 共␦r兲ei. Then, we have
036107-2
PHYSICAL REVIEW E 77, 036107 共2008兲
SYNCHRONIZATION IN NETWORKS OF NETWORKS: THE ...
z = g共s兲 兺 k⬘e−i␣⬘z⬘ . ⬘
Now define the complex quantity ¯k⬘ = k⬘e−i␣⬘. Using the Kronecker delta ␦⬘, the above equation can be written
兺 ⬘
冉
¯k − ⬘
␦⬘
g共s兲
冊
tions; we set ⌬ = ⌬ = ⌬ and ⍀ = ⍀ = ⍀. 共A more generic example will follow.兲 Denoting D = det共K兲 and T = tr共K兲, we separate the real and imaginary parts of Eq. 共7兲 to obtain D2 − 2⌬T + 4⌬2 − 4共v + ⍀兲2 = 0, 共v + ⍀兲共4⌬ − T兲 = 0.
z⬘ = 0.
One solution to these equations is
In matrix notation for the case of two populations labeled and , this is
冢
¯k −
1 g共s兲
¯k
¯k ¯k − 1 g共s兲
冣
冉冊
z = 0. z
共5兲
This equation has a nontrivial solution if the determinant of the matrix is zero. This condition determines the growth rate s in terms of , ¯k⬘, and the parameters of the distributions G共兲 关via g共s兲兴. To illustrate the resulting behavior, we note that g共s兲 can be evaluated in closed form for a Cauchy-Lorentz distribution G 共 兲 =
1 ⌬ , 共 − ⍀兲2 + ⌬2
共6兲
where ⍀, the mean frequency of population , and the half width at half maximum ⌬ are both real, and ⌬ is positive. One obtains
ⴱ = ⌬
vⴱ = − ⍀,
Using this, the determinant condition for the two-population case reduces to ¯ − 2共s + ⌬ + i⍀ 兲兴关k ¯ − 2共s + ⌬ + i⍀ 兲兴 − ¯k ¯k = 0. 关k For simplicity, we set the phase angles ␣⬘ to zero for the remainder of this paper 关14兴. In this case, the matrix elements ¯k ¯ ⬘ are purely real, so that k⬘ = k⬘. The determinant condition then becomes 关k − 2共s + ⌬ + i⍀兲兴关k − 2共s + ⌬ + i⍀兲兴 共7兲
Notice that if = 0, then s = −⌬ − i⍀, indicating that the incoherent state is stable for zero coupling 共since −⌬ is negative兲. We imagine increasing 共or decreasing兲 until s crosses the imaginary axis at a critical value ⴱ. At this point, the incoherent state loses stability and coherent collective behavior emerges in the ensemble. Thus the critical value ⴱ can be determined from the determinant condition by setting s = iv, where v is real 共so that the perturbations are marginally stable兲, and equating the real and imaginary parts of the left side of Eq. 共7兲 to zero. This results in two equations which can be solved simultaneously for the two 共real兲 unknowns and v. For our first example, we choose two identical popula-
冉
冊
T ⫾ 冑T2 − 4D , D
which is valid for T2 ⱖ 4D, since is assumed to be real. Notice that the appropriate solution as D → 0 共using the negative sign for T ⬎ 0 and the positive sign for T ⬍ 0兲 is vⴱ = − ⍀,
ⴱ =
2⌬ , T
as can be verified by setting D = 0 in Eq. 共8兲 directly. Another solution is vⴱ = − ⍀ ⫾
⌬ 冑4D − T2, T
ⴱ =
4⌬ . T
Finally, setting T = 0 in Eq. 共8兲 yields vⴱ = −⍀ and ⴱ ⌬ = ⫾ 冑2−D for D ⬍ 0, and no solution for D ⱖ 0. These results are summarized in Table I. Thus the critical values ⴱ are determined by T, D, and ⌬. To illustrate this result, we begin by discussing a particular example. Consider the matrix K=
1 = 2共s + ⌬ + i⍀兲. g 共 兲
− 2kk = 0.
共8兲
冉 冊 1 −1 1
0
,
which has trace T = 1 and determinant D = 1. This corresponds to case 2 in Table I, from which we find that ⴱ = 4⌬ / T. Figure 1 shows the results of a numerical simulation of two populations of 10 000 oscillators each, using ⌬ = 1. The order parameters r vs are shown, and we can see that the oscillator populations are incoherent for values below the predicted critical value ⴱ = 4, and that they grow increasingly synchronized as is increased beyond ⴱ. Next, we examined eight different connectivity matrices K that were chosen to sample the various regions in T − D space. Table II shows these matrices and the corresponding value共s兲 of ⴱ. These predictions were tested by numerically calculating the order parameters r for a range of coupling TABLE I. Solutions to Eq. 共8兲 for two identical populations. D = det共K兲 and T = tr共K兲, where K is the connectivity matrix 关Eq. 共1兲兴, and ⌬ is the width parameter in Eq. 共6兲. Case
Condition
vⴱ
1
T2 ⬎ 4D
−⍀
2
T ⱕ 4D
3
D=0
036107-3
2
−⍀ ⫾ T 冑4D − T −⍀ ⌬
4
T = 0, D ⬍ 0
−⍀
5
T = 0, D ⱖ 0
no solution
ⴱ ⌬共 2
T⫾冑T2−4D
D
兲
4⌬
T
2⌬
T
2⌬
⫾ 冑−D no solution
PHYSICAL REVIEW E 77, 036107 共2008兲
BARRETO et al.
TABLE II. Connectivity matrices K chosen to sample T-D space. Case A
B
C
D FIG. 1. Numerical calculation of the order parameters 共쎲 , △兲 vs for case A in Table II. The vertical line corresponds to the predicted value ⴱ = 4. The data point nearest ⴱ is at = 4.15.
values 关15兴. Results are shown in Fig. 2. The various cases are located by letter in the T − D plane according to the trace and determinant of the matrices, and the corresponding inset shows the numerically calculated order parameters plotted vs , with the predicted critical coupling indicated by a vertical line at = ⴱ. For example, the inset to case A 关with 共T , D兲 = 共1 , 1兲兴 corresponds to Fig. 1, which was discussed above. It can be seen that for all cases, the onset of synchronization occurs as predicted. Note that more than one prediction for ⴱ may be specified by our analysis 共see Table I兲. The solutions closest to = 0 are the relevant ones, because we expect the incoherent state to lose stability once the first ⴱ solution is encountered. There are two possible cases depending on the sign of D. First, if the two solutions have the same sign, then there is only one critical value ⴱ 共which may be positive or negative depending on the sign of the trace兲 for the onset of synchronization. This occurs for D ⬎ 0 and T ⫽ 0, as can be seen in
E
F
G
H
T
D
ⴱ
1 −1 1 0
1
1
4
−2 −3 1 1
−1
1
−4
3 1 −3.5 −1
2
0.5
2共2 − 冑2兲 = 1.172
−3 1 −3.5 1
−2
0.5
2共−2 + 冑2兲 = −1.172
−1 −1 1 2
1
−1
−共1 ⫾ 冑5兲 = −3.236, 1.236
1 1 −1 −2
−1
−1
1 ⫾ 冑5 = −1.236, 3.236
2 1 −3 −2
0
−1
⫾2
1 −1 2 −1
0
1
none
Matrix
共 兲 共 兲
共 兲 共 兲 共 兲 共 兲 共 兲 共 兲
Fig. 2. 共Interestingly, for D ⬎ 0 and T = 0, synchronization does not occur for any .兲 The other case occurs for D ⬍ 0, for which the two ⴱ solutions have opposite signs. In this case, there are two critical values ⴱ for the onset of synchronization—one on either side of = 0. This can also be observed in Fig. 2. In the more general case in which the various populations have different natural frequency distributions, it is not typi-
FIG. 2. Numerical simulations using the matrices listed in Table II for identical populations. The letters indicate the placement of each case in the T-D plane, and the corresponding insets show numerical calculations of the order parameters 共쎲 , △兲 vs for that case 共in all cases, = 0 is in the center of the horizontal axis兲. Vertical lines in the insets indicate the predicted value共s兲 ⴱ listed in Table II for the onset of coherent collective behavior. In all cases, we used ⌬ = 1. Note that for D ⬎ 0, there is one value of ⴱ, whose sign corresponds to the sign of the trace T. If D ⱖ 0 and T = 0, then synchronization is not possible for any coupling strength. For D ⬍ 0, there are two values of ⴱ: one positive, and one negative.
036107-4
PHYSICAL REVIEW E 77, 036107 共2008兲
SYNCHRONIZATION IN NETWORKS OF NETWORKS: THE ... 6
6
4
Im(η)
Im(η)
4 2 0 −2
−4 −4
−3
v
−2
−1
−5
−4
−6
−5
−4
v
−3
−2
−1
−3
−2
−1
6
0
4
−1
Re(η)
Re(η)
−6
(a)
1
−2
2 0
−3 −4
0 −2
(a)
(b)
2
−4
−3
v
−2
−2
−1
(b)
FIG. 3. Case E, different populations. The upper panel shows Im共共1,2兲兲, with the vertical lines identifying roots at v1 = −4.024 and v2 = −1.722. The lower panel shows Re共共1,2兲兲; values at the roots found above are indicated by horizontal lines, yielding 共1兲 ⴱ = 0.515 and 共2兲 ⴱ = −2.809. Thus we expect synchronization to occur at these values as is either increased or decreased away from zero.
cally possible to describe the onset of synchronization in terms of the determinant and trace of the coupling matrix K = k⬘ alone. We now consider this situation, but retain the Cauchy-Lorentz form of the natural frequency distributions for convenience. We manipulate Eq. 共7兲 as follows. Let s = iv 共i.e., purely imaginary, to consider the marginally stable case兲 and define a = ⌬ + i共v + ⍀兲 and b = ⌬ + i共v + ⍀兲. We obtain
v
FIG. 5. Case A, different populations. Panels as in Fig. 3. Note that the standard branch cut leads to discontinuities and the occurrence of two roots for Im共共1兲兲共v兲 共solid lines, upper panel兲, and none for Im共共2兲兲共v兲 共dotted lines, upper panel兲. From the lower panel we find ⴱ = 2.189, 4.501, taking care to obtain these from Re共共1兲兲 共solid line, lower panel兲. We expect to find synchronization onset at the smaller of these values.
There are two unknowns in this equation: and v. Equation 共9兲 is a quadratic equation in with complex coefficients,
and we can easily obtain two complex solutions 共1,2兲 as functions of v. Since the critical values 共1,2兲 must be real, ⴱ we solve for the roots of Im共共1,2兲兲, and evaluate Re共共1,2兲兲 at . these roots. This yields the possible critical values 共1,2兲 ⴱ Typically, these steps must be performed symbolically and/or numerically; we used MATLAB® 关16兴. As before, the values of 共1,2兲 that are closest to zero 共on either side兲 are the relⴱ evant values. To illustrate, we choose two populations with CauchyLorentz natural frequency distributions 关Eq. 共6兲兴 with parameters ⌬ = 1, ⍀ = 2, ⌬ = 0.5, and ⍀ = 4. We consider the
FIG. 4. Case E, different populations. Calculations of the order parameters vs confirm that the onset of synchronization occurs at ⴱ = −2.809 and 0.515 共vertical lines兲, as predicted in Fig. 3.
FIG. 6. Case A, different populations: The onset of synchronization occurs at ⴱ = 2.189 共vertical line兲, as predicted in Fig. 5. No synchronization is observed for smaller values of 共not shown兲.
D2 − 2共bk + ak兲 + 4ab = 0.
共9兲
036107-5
PHYSICAL REVIEW E 77, 036107 共2008兲
6
3
4
2
Im(η)
Im(η)
BARRETO et al.
2 0
−6 (a)
1 0 −1
−5
−4
−3
v
−2
−1
0
−2
1 (a)
−4
−3
−4
−3
−2
−1
0
−2
−1
0
v
6
3
4
2
2
Re(η)
Re(η)
8
0 −2 −4
(b)
−6
FIG. 7. Case = −1.429, 5.000.
−5
H,
−4
−3
v
different
−2
−1
0
populations.
0 −1
1
We
1
find
ⴱ
same K matrices as before, i.e., those listed in the second column of Table II. Case E is straightforward; the analysis is illustrated in Fig. 3 and the numerical verification is in Fig. 4. Note that since Eq. 共9兲 has complex coefficients, obtaining typically involves taking the square root of a complex number. Therefore one must be mindful of branch cuts when obtaining symbolic and/or numerical solutions. This is important in the analysis for case A, shown in Figs. 5 and 6. Finally, case H, which exhibits no synchronization for identical populations for any value of , does show synchronization in the present case. The analysis is shown in Figs. 7 and 8. We close by giving an example with three different populations. We choose the same Cauchy-Lorentz distributions as above, and add a third with ⌬ = 1 / 3 and ⍀ = 1. We use the following K matrix:
FIG. 8. Case H, different populations. Synchronization occurs at
ⴱ = −1.429 and 5.000 共vertical lines兲, as predicted in Fig. 7.
(b)
v
FIG. 9. 共Color online兲 Three populations. The imaginary and real parts of 共1,2,3兲 are plotted to illustrate the discontinuities due to branch cuts. The analysis proceeds as in the previous cases. Because of branch cuts, two roots occur on 共1兲, one on 共2兲, and none on 共3兲. Evaluating the real parts appropriately, we find ⴱ = −0.891, −0.564, 2.303.
冢
−1
1
1
1
−1
1
1
1
−1
冣
.
The procedure for deriving ⴱ proceeds as above, except that Eq. 共7兲 is replaced by a third-degree polynomial in . Figure 9 shows the imaginary and real parts of the three solutions. 共Note that the branch cuts are more complicated.兲 The predicted onset of synchronization was verified, as shown in Fig. 10.
FIG. 10. Three populations. Synchronization occurs at ⴱ = −0.564 and 2.303, as predicted in Fig. 9.
036107-6
PHYSICAL REVIEW E 77, 036107 共2008兲
SYNCHRONIZATION IN NETWORKS OF NETWORKS: THE ...
ACKNOWLEDGMENTS
In conclusion, we have described how to determine the onset of coherent collective behavior in systems of interacting Kuramoto systems, i.e., systems of interacting populations of phase oscillators with both node and coupling heterogeneity.
E.B. was supported by NIH Grant No. R01-MH79502; E.O. was supported by ONR 共Physics兲 and by NSF Grant No. PHY045624.
关1兴 R. Milo et al., Science 298, 824 共2002兲. 关2兴 M. E. J. Newman and M. Girvan, Phys. Rev. E 69, 026113 共2004兲. 关3兴 A. Arenas, A. Díaz-Guilera, and C. J. Pérez-Vicente, Phys. Rev. Lett. 96, 114102 共2006兲; Physica D 224, 27 共2006兲. 关4兴 M. Kurant and P. Thiran, Phys. Rev. Lett. 96, 138701 共2006兲. 关5兴 X. Wang, L. Huang, Y.-C. Lai, and C. H. Lai, Phys. Rev. E 76, 056113 共2007兲. 关6兴 C. Zhou, L. Zemanová, G. Zamora, C. C. Hilgetag, and J. Kurths, Phys. Rev. Lett. 97, 238103 共2006兲. 关7兴 M. Girvan and M. E. J. Newman, Proc. Natl. Acad. Sci. U.S.A. 99, 7821 共2002兲. 关8兴 R. Milo et al., Science 303, 1538 共2004兲. 关9兴 Y. Kuramoto, in International Symposium on Mathematical Problems in Theoretical Physics, edited by H. Araki, Lecture Notes in Physics Vol. 39 共Springer, New York, 1975兲, 420; Chemical Oscillators, Waves and Turbulence 共Springer, Berlin, 1984兲. 关10兴 A. T. Winfree, The Geometry of Biological Time 共Springer, New York, 1980兲; S. H. Strogatz, Physica D 143, 1 共2000兲; J. A. Acebrón et al.; Rev. Mod. Phys. 77, 137 共2005兲; K. Wiesenfeld and J. W. Swift, Phys. Rev. E 51, 1020 共1995兲; I. Z. Kiss, Y. Zhai, and J. L. Hudson, Science 296, 1676 共2002兲. 关11兴 A similar system of two asymmetrically interacting populations with a particular form of coupling matrix K was consid-
ered in E. Montbrio, J. Kurths, and B. Blasius, Phys. Rev. E 70, 056125 共2004兲. The formulation in the current paper is more general in two important aspects: K is arbitrary, and we allow for any number of interacting populations. Our system is similar to that studied in J. G. Restrepo, E. Ott, and B. R. Hunt, Chaos 16, 015107 共2006兲. However, our formulation permits different natural frequency distributions for each population. The bracketed expression is valid for Re共s兲 ⬎ 0 and may be analytically continued into the region where Re共s兲 ⱕ 0. See E. Ott, P. So, E. Barreto, and T. Antonsen, Physica D 173, 29 共2002兲; E. Ott, Chaos in Dynamical Systems, 2nd ed. 共Cambridge University Press, Cambridge, England, 2002兲, p. 240. Many interesting states require ␣⬘ ⫽ 0, such as the chimera state observed in D. M. Abrams and S.H. Strogatz, Int. J. Bifurcation Chaos Appl. Sci. Eng. 16, 21 共2006兲. Simulations used fourth-order Runge-Kutta with a time step of 0.01 s, N = 10 000 or 50 000, and parameters as noted. The system was initialized in the incoherent state and an initial transient was discarded. The order parameters were then averaged over the subsequent 10 s. Because the standard deviations over this interval were small, no error bars were plotted. MATLAB®, The MathWorks, Inc., Natick, MA, http:// www.mathworks.com/
关12兴
关13兴
关14兴
关15兴
关16兴
036107-7
