Nonstationary pattern in unsynchronizable complex networks
Nonstationary pattern in unsynchronizable complex networks
Xingang Wang,1, 2 Meng Zhan,3 Shuguang Guan,1, 2 and Choy Heng Lai4, 2
1Temasek Laboratories, National University of Singapore, 117508 Singapore
2Beijing-Hong Kong-Singapore Joint Centre for Nonlinear & Complex Systems (Singapore),
National University of Singapore, Kent Ridge, 119260 Singapore
3Wuhan Institute of Physics and Mathematics, Chinese Academy of Sciences, Wuhan 430071, China
4Department of Physics, National University of Singapore, 117542 Singapore
(Dated: November 12, 2018)
Pattern formation and evolution in unsynchronizable complex networks are investigated. Due to the asym-
metric topology, the synchronous patterns formed in complex networks are irregular and nonstationary. For
coupling strength immediately out of the synchronizable region, the typical phenomenon is the on-off inter-
mittency of the system dynamics. The patterns appeared in this process are signatured by the coexistence of a
giant cluster, which comprises most of the nodes, and a few number of small clusters. The pattern evolution
is characterized by the giant cluster irregularly absorbs or emits the small clusters. As the coupling strength
leaves away from the synchronization bifurcation point, the giant cluster is gradually dissolved into a number of
small clusters, and the system dynamics is characterized by the integration and separation of the small clusters.
Dynamical mechanisms and statistical properties of the nonstationary pattern evolution are analyzed and con-
ducted, and some scalings are newly revealed. Remarkably, it is found that the few active nodes, which escape
from the giant cluster with a high frequency, are independent of the coupling strength while are sensitive to the
bifurcation types. We hope our findings about nonstationary pattern could give additional understandings to
the dynamics of complex systems and have implications to some real problems where systems maintain their
normal functions only in the unsynchronizable state.
PACS numbers: 89.75.-k, 05.45.Xt
I. INTRODUCTION
Synchronization of complex networks has aroused many in-
terest in nonlinear science since the discoveries of the small-
world and scale-free properties in many real and man-made
systems [1, 2]. In this study, one important issue is to ex-
plore the inter-dependent relationship between the collective
behaviors of the complex systems and their underlying topolo-
gies. In particular, many efforts have been paid to the con-
struction of optimal networks, and a number of factors which
have important affections to the synchronizability of complex
networks have been gradually disclosed. Now it is known that
random networks, due to their small average distances, are
generally more synchronizable than regular networks [3, 4];
and scale-free networks, with weighted and asymmetric cou-
plings, can be more synchronizable than homogeneous net-
works [5, 6]. In these studies, the standard method employed
for synchronization analysis is the master stability function
(MSF), where the network synchronizability is estimated by
an eigenratio calculated from the coupling matrix, and system
which has a smaller eigenratio is believed to be more synchro-
nizable than that of larger eigenratio [7]. Inspired by this, to
improve the network synchronizability, the only task seems to
be upgrading the coupling matrix so as to decrease the eigen-
ratio, either by changing the network topology [4] or by ad-
justing the coupling scheme [5, 6].
The MSF method, while bringing great convenience to the
analysis, overlooks the temporal, local property of the sys-
tem and reflects only partial information about the system
dynamics. Specifically, from MSF we only know ultimately
whether the network is globally synchronizable or unsynchro-
nizable, but do not know how the global synchronization is
reached if the network is synchronizable, or what’s the pat-
tern and how it evolves if the network is unsynchronizable.
These evolution details, or the transient behavior in system
development, contain rich information about the system dy-
namics and may give additional insights to the organization
of complex systems. For instance, the recent studies about
synchronization transition have shown that, in the unsynchro-
nizable states, heterogeneous networks are more synchroniz-
able (have a higher degree of coherence) than homogeneous
networks at small couplings, while at larger couplings the op-
posite happens [8]. This crossover phenomenon of network
synchronizability are difficult to understood if we only look
at the final state of the system, but are straightforward if we
look at the transient behaviors of their evolutions [8]. Besides
revealing the synchronization mechanisms, the transient be-
havior of network synchronization can also be used to detect
the topological scales and hierarchical structures in the real
systems, e.g., the detection of cluster structures in social and
biological networks [9, 10]. However, despite of its theoreti-
cal and practical significance, the study of transient dynamics
of complex networks is still at its infancy and many questions
remain open, say for example, the pattern evolution of unsyn-
chronizable complex networks.
Pattern formation in unsynchronizable but near-
synchronization networks has been an important issue
in studying the collective behavior of regular networks
[11, 12]. By setting the coupling strength nearby the
synchronization bifurcation point, the system state shares
both the dynamical properties of the synchronizable and
unsynchronizable states: a state of high coherence but is not
synchronized. The bifacial dynamical property makes this
state a natural choice in investigating the transition process of
networks synchronization. Previously studies about regular
http://arxiv.org/abs/0704.0892v1
networks, say for example the lattices [11], have shown that,
when the coupling strength is slightly out of the synchroniz-
able region, although global synchronization is unreachable,
nodes are still synchronized in a partial sense. That is, nodes
are self-organized into a number of synchronous clusters.
The distribution of these clusters, also called the synchronous
pattern, is determined by a set of factors such as the coupling
strength, the system size and the coupling scheme. As the
coupling strength leaves away from the bifurcation point, the
pattern structure becomes more and more complicated and the
system coherence will be decreased, and finally reaches the
turbulence state. It is worthy of note that the patterns arisen
in regular networks have two common properties: spatially
symmetric and temporally stationary. More specifically, the
contents of each cluster are fixed and the clusters are of
translation symmetry in space. For this reason, we say that
the patterns formed in regular networks are symmetric and
stationary. These two properties, as have been discussed in
the previous studies [11, 12], are rooted in the symmetric
topology of the regular networks. This makes it interesting
to ask the following question: how about the patterns in
unsynchronizable complex networks?
Different to the regular networks, in complex networks we
are not able to find any symmetry from their topologies. The
asymmetric topology, according to the pattern analysis devel-
oped in studying regular networks [11], will induce two sig-
nificant changes to the patterns: 1) the synchronous clusters, if
they exist, will be asymmetric; and 2) all the possible patterns,
including the one of global synchronization, are linearly un-
stable under small perturbations. In other words, the patterns
formed in complex networks are expected to be asymmetric
and nonstationary. Our mission of this paper is just to under-
stand and characterize the nonstationary patterns arisen in the
development of complex networks. Specifically, we are trying
to investigate the following questions: 1) is there any pattern
arises during the system evolution? 2) the pattern is stationary
or nonstationary? if nonstationary, how is it evolving and how
is it reflected from the system dynamics? 3) What happens
to the pattern properties during the transition of network syn-
chronization? and 4) How the coupling strength and bifurca-
tion type affect the pattern properties? By investigating these
dynamical and statistical properties, we wish to have a global
understanding to the dynamics of unsynchronizable complex
networks.
Our main findings are: 1) for coupling strength immedi-
ately outside of the synchronizable region, the system dynam-
ics undergoes the process of on-off intermittency. That is,
most of the time the system stays on the global synchroniza-
tion state (the ”off” state) but, once in a while, it develops
into a breaking state (the ”on” state) which is composed by
a giant cluster and a few number of small clusters (hereafter
we call it the giant-cluster state). As the system develops,
the giant cluster changes its shape by absorbing or emitting
the small clusters, leading to the ”off” or ”on” states, respec-
tively; 2) the few active nodes which escape from the giant
cluster with the high frequencies are coupling-strength inde-
pendent but are bifurcation-type dependent. That is, in the
neighboring region of a fixed bifurcation point, the locations
of these active nodes do not change with the coupling strength;
if we change the coupling strength from nearby another bi-
furcation point (the two bifurcation points will be explained
later), their locations will be totally changed; 3) as coupling
strength leaves away from the bifurcation point, the giant clus-
ter is gradually dissolved and more small clusters are gener-
ated from it. Eventually, the giant cluster disappears and the
pattern is composed by only the small clusters (hereafter we
call it the scattering–cluster state). During the course of sys-
tem evolution, each small cluster may either increase its size
by integrating with other small clusters or decrease its size
by breaking to even small clusters, but it can never reach to
the global synchronization state; 4) besides the giant cluster,
the giant- and scattering-cluster states are also distinct in their
small clusters. For giant-cluster state the size of the small
clusters follows a power-law distribution, while for scattering-
cluster state it follows a Gaussian distribution.
The rest of the paper is going to be arranged as follows. In
Sec. II we will give our model of coupled map network and,
based on the method of MSF, point out the two bifurcation
points and the transition areas that we are going to study with.
In Sec. III we will employ the method of finite-time Lyapunov
exponent to predict and describe the intermittent system dy-
namics in the bifurcation regions. Direct simulations about
on-off intermittency will be presented in Sec. IV. By intro-
ducing the method of temporal phase synchronization, in Sec.
V we will investigate in detail the dynamical and statistical
properties of the nonstationary pattern. Meanwhile, proper-
ties of the giant- and scattering states will be compared and
the transition between the two states will be conducted. In
Sec. VI we will discuss the phenomenon of active nodes and
investigate their dependence to the network properties. Dis-
cussions and conclusions about pattern evolution in complex
networks will be presented in Sec. VII.
II. COUPLED MAP NETWORKS AND THE TWO
BIFURCATION POINTS
Our model of coupled map network is of the following form
xi(t+ 1) = F(xi(t))− ε
Gi,jH [f(xj(t))] . (1)
where xi(t + 1) = F(xi(t)) is a d-dimensional map repre-
senting the local dynamics on node i, ε is a global coupling
parameter, G is Laplacian matrix representing the couplings,
and H is a coupling function. To facilitate our analysis, we
adopt the following coupling scheme [13]:
Gi,j = −
Ai,jk
j=1 Ai,jk
, (2)
for j 6= i and Gi,i = 1, with ki the degree of node i and
A the adjacent matrix of the network: Ai,j = 1 if node i and
j are connected and Ai,j = 0 otherwise. In comparison with
the traditional coupling schemes, one important advantage we
benefit from this coupling scheme is that the synchronizabil-
ity of the network, i.e. the eigenratio of the coupling matrix
described in Eq. [2], can be easily adjusted by the parameter
β, while the network topology is kept unchanged. This advan-
tage brings many convenience in network selection since for
a given network topology, even though it is unsynchronizable
under the traditional schemes, can now be synchronizable by
adjusting β in Eq. [2]. This convenience is of particular im-
portance when our studies of network dynamics are focused
on the bifurcation regions, where the network synchronizabil-
ity should be deliberately arranged in order to demonstrate
both the two types of bifurcations. We note that the adop-
tion of Eq. [2] is only for the purpose of convenient analy-
sis, the findings we are going to report are general and can
also be observed by other coupling schemes given the net-
work is properly prepared. In practice, we use logistic map
F(x) = 4x(1−x) as the local dynamics and adopt H(x) = x
as the coupling function.
We first locate the two bifurcation points of global syn-
chronization. The linear stability of the global synchroniza-
tion state {xi(t) = s(t), ∀i} is determined by the correspond-
ing variational equations, which can be diagonalized into N
blocks of form
y(t+ 1) = [DF(s) + σDH(s)] y(t), (3)
with DF(s) and DH(s) the Jacobian matrices of the corre-
sponding vector functions evaluated at s(t), and y represents
the different modes that are transverse to the synchronous
manifold s(t). We have σ(i) = ελi for the ith block, i =
1, 2, ..., N , and λ1 = 0 ≤ λ2 ≤ ... ≤ λN are the eigenvalues
of matrix G. The largest Lyapunov exponent Λ(σ) of Eq. [3],
known as the master stability function (MSF) [7], determines
the linear stability of the synchronous manifold s(t). In par-
ticular, the synchronous manifold is stable if only Λ(ελi) < 0
for each i = 2, ..., N . The set of Lyapunov exponents Λ(ελi)
govern the stability of the synchronous manifold in the trans-
verse spaces, and a positive value of Λ(ελi) represents the
loss of the stability in the transverse space of mode i. It was
found that for a large class of chaotic systems, Λ(σ) < 0
is only fulfilled within a limit range in the parameter space
σ ∈ (σ1, σ2). This indicates that, to make the global synchro-
nization state linearly stable, all the eigenvalues λi should be
contained within range (σ1, σ2), i.e., λN/λ2 < σ2/σ1. For
the logistic map employed here, it is not difficult to prove that
σ1 = 0.5 and σ2 = 1.5. Therefore, to achieve global syn-
chronization, the coupling matrix G should be designed with
eigenratio R ≡ λN/λ2 < σ2/σ1 = 3 = Rc.
Besides the condition of R < Rc, to guarantee the synchro-
nization, we also need to set the coupling strength in a proper
way: either small or large couplings may deteriorate the syn-
chronization. If ε < ε1 = σ1/λ2, the couplings are too weak
to restrict the node trajectories to the synchronous manifold;
while if ε > ε2 = σ2/λN , the couplings will be too strong
and actually act as large perturbations to the synchronization
manifold. Therefore, to achieve the global synchronization,
we also require ε1 < ε < ε2. The two critical couplings
ε1 and ε2, which are named as the long-wave (LW) [14] and
short-wave (SW) bifurcations [15] respectively in the studies
of regular networks, thus stand as the boundaries of the syn-
chronizable region. Our studies about network synchroniza-
0.80 0.85 0.90 0.95 1.00
 = 0.95127
 = 0.83488
FIG. 1: For scale-free network of N = 1000 nodes and of average
degree < k >= 8, a schematic plot on the generation of the two
bifurcation points as a function of the coupling strength. The long-
wave bifurcation occurs at about ε1 ≈ 0.83488 which is determined
by the condition ελ2 = σ1 = 0.5 (the lower line). The short-wave
bifurcation occurs at about ε2 ≈ 0.95127 which is determined by
the condition ελN = σ2 = 1.5 (the upper line).
tion will be focused on the neighboring regions of the two
bifurcation points, i.e., the region of ε . ε1 or ε & ε2.
By the standard BA growth model [1], we construct a scale-
free network of 103 nodes and of average degree 〈k〉 = 8. By
setting β = 2.5 in Eq. [2], we have λ2 ≈ 0.6 and λN ≈ 1.58.
Because of R = λN/λ2 ≈ 2.6 < Rc, the network is glob-
ally synchronizable. Also, because of λ2 > σ1 and λN > σ2,
both the two bifurcations can be realized by adjusting the cou-
pling strength within range ε ∈ (0, 1). In specific, when
ε < ε1 ≈ 0.835, we have ελ2 < σ1 and ελN < σ2, the
synchronous manifold loses its stability at the lower boundary
of the synchronizable region and LW bifurcation occurs; and
when ε > ε2 ≈ 0.95, we have ελ2 > σ1 and ελN > σ2, the
synchronous manifold loses its stability at the upper bound-
ary of the synchronizable region and SW bifurcation occurs
[Fig. 1]. In the following we will fix the network topology
and the parameter β, while generating the various patterns by
changing the coupling strength ε nearby the two bifurcation
points.
III. FINITE-TIME LYAPUNOV EXPONENT
Before direct simulations, we first give a qualitative de-
scription (prediction) on the possible system dynamics in bi-
furcation regions. To concrete our analysis, in the following
we will only discuss the situation of SW bifurcation (ε . ε1),
while noting that the same phenomena can be found at the
LW bifurcation as well (ε & ε2) . In preparing the unsynchro-
nizable states, we only let Λ(λ2) be slightly puncturing into
the unstable region, while keeping all the other exponents still
staying in the stable region, i.e., Λ(λ2) & 0 and Λ(λi) < 0
for i = 3, ...N . With this setting, the synchronous manifold
is only desynchronized in the transverse space of mode 2. As
such, the system possesses only two positive Lyapunov expo-
nents, one is Λ(λ0) which is associated to the synchronous
manifold itself and another one is Λ(λ2). Noticing that Λ(λ)
are asymptotic averages, and, as so, they account only for the
global stability properties, but do not warrant the possible co-
herent sets arising in the system evolutions. These coherent
sets, for regular networks, refer to the stationary, symmetric
patterns to which the system finally develops. While for com-
plex networks, these sets can be the temporal, irregular clus-
ters emerged in the process of system evolution.
In the region of ε . ε1, although global synchronization is
unreachable, the system may still keep with the high coher-
ence due to the existence of the synchronous clusters. Espe-
cially, there could be some moments at which all the trajecto-
ries are restrained to a small region in the phase space, very
close to the situation of global synchronization. This vary-
ing system coherence, however, can not be reflected from the
asymptotic value Λ(λ). To characterize the variation, we need
to employ some new quantities which are able to capture the
temporal behavior of system. One of such quantities is the
finite-time Lyapunov exponent (FLE), a technique developed
in studying chaos transition in nonlinear science [16]. In stead
of asymptotic average, FLE measures the diverging rate of
nearby trajectories only in a finite time interval T .
t=(i−1)T
lnDH(s(t)). (4)
As our studies are focused on the situation of one-mode
desynchronization, the stability of the synchronous manifold
and the temporal behavior that it displays are therefore ex-
pected to be more reflected from the variation of Λ2,i, the
FLE that associates with mode 2. With ε = 0.83 and
T = 100, we plot in Fig. 2 the time evolution of Λ2,i. It
is found that, although with a positive asymptotic value about
〈Λ2,i〉 ≈ 6 × 10
−3, the instant value of Λ2,i penetrates to the
negative region at a high frequency. According to the different
signs of Λ2,i, the system evolution is divided into two types of
intervals: the divergent interval and the contractive interval.
In the divergent intervals we have Λ2,i > 0 and the system
dynamics is temporarily dominated by the divergence of the
node trajectories from the synchronous manifold; while in the
contractive intervals we have Λ2,i < 0 and the system dynam-
ics is temporarily dominated by the convergence of the node
trajectories to the synchronous manifold.
The variation of Λ2,i, reflected on the process of pattern
evolution, characterizes the travelling property of the sys-
tem dynamics among the neighboring regions of two different
kinds of states: the desynchronization state and the synchro-
nization state. In Fig. 2, the minimum value of Λ2,i is about
−0.07, during this contractive interval the node trajectories
will converge to the synchronous manifold by an amount of
eminΛ2,iT ≈ e−7 ≈ 10−3 on average. Assuming that before
entering this interval the average distance between the node
trajectories is ∆ (for logistic map we always have ∆ < 1),
then at the end of this interval the average distance is de-
creased to ∆ × 10−3, a small value which is usually over-
shadowed by noise in practice. Due to this small distance,
0 400 800 1200 1600 2000
-0.08
-0.04
FIG. 2: For ε = 0.83 in Fig. 1, the time evolution of the finite-
time Lyapunov exponent Λ2,i calculated on intervals of length T =
100. It is observed that, while having the positive asymptotic value
〈Λ2,i〉 > 0, the temporal value of Λ2,i is penetrating into the negative
region frequently.
the system can be practically regarded as already reached the
synchronization state. On the other hand, if the system enters
a divergent interval, the node trajectories will diverge from
each other and, at the end of this interval, their average dis-
tance will be increased by an order of 103. This large distance
will deteriorate the ordered trajectories (or the high coherence
of the system dynamics) that achieved during the contractive
intervals, and leading to the incoherent, breaking state. The
pattern of the breaking state, however, is not unique. Depend-
ing on the initial conditions and the divergence intervals, the
pattern may assume the different configurations. Therefore,
based on the observations of Λ2,i [Fig. 2] the dynamics of
unsynchronizable networks can be intuitively understood as
an intermittent travelling among the synchronization state and
the different kinds of desynchronization states.
IV. ON-OFF INTERMITTENCY DESCRIBED BY
COMPLETE SYNCHRONIZATION
We now investigate the system dynamics by direct simula-
tions. To implement, we first prepare the system to be staying
on the synchronization state. This can be achieved by adopt-
ing a large coupling strength from the synchronizable region,
i.e. ε1 < ε < ε2. After synchronization is achieved, we then
decrease ε to a value slightly below the bifurcation point ε1
and, in the meantime, an instant small perturbation is added
on each node. In practice, we take i.i.d (independent iden-
tically distributed) noise of strength 1 × 10−5 as the pertur-
bations. After this, we release the system and let it develop
according to Eq. (1). Since ε < ε1, the synchronization state
is unstable and, triggered by the noises, the node trajectories
begin to diverge from each other. The divergent trajectories,
however, will frequently visit the neighborhood of the syn-
chronous manifold, especially during those contractive inter-
vals of small Λ2,i [Fig. 2]. The intermittent system dynamics
is plotted in Fig. 3(a), where the average trajectory distance
∆X = 1
i=1 xi − ~x is plotted as a function of time. As
we have predicted from LLE, the system dynamics indeed un-
dergoes an intermittent process. To characterize the intermit-
tency, we plot in Fig. 3(b) the laminar-phase distribution of
the ∆X sequence plotted in Fig. 3(a). It is found that the lam-
inar length τ (the time interval between two adjacent bursts
of amplitude ∆X(t) > 10−3) and the probability p(τ) for it
to appear follow a power-law scaling p(τ) ∼ τ−γ . The fitted
exponent is about γ ≈ −1.5 ± 0.05, with a fat tail at large τ
due to the finite simulating time.
In chaos theory, intermittent process of laminar-phase expo-
nent −3/2 is classified as the ”on-off” intermittency, a typical
phenomenon observed in dynamical systems with a symmet-
ric invariant set [17]. On-off intermittency is also reported
in chaos synchronization of regular networks, where the in-
variant set refers to the synchronous manifold, and the ”off”
state refers to the long stretches that the system dynamics is
staying nearby the synchronous manifold and the ”on” state
refers to the short bursts that the system dynamics is staying
away from the synchronous manifold. Therefore, in terms of
laminar-phase distribution, the intermittency we have found
in complex networks [Fig. 3] has no difference to the that of
the regular networks, despite of the drastic difference between
their topologies. We have also investigated the transition be-
havior of the averaged distance 〈∆X(t)〉 nearby the bifurca-
tion points. As shown in Fig. 3(c), a linear relation between
〈∆X(t)〉 and ε is found in the region of ε . ε1. This linear
transition of the system performance, again, is consistent with
the transition of regular networks [18]. Therefore, in terms of
complete synchronization, the on-off intermittency we have
found in complex networks has no difference to that of the
regular networks.
V. PATTERN EVOLUTION IN COMPLEX NETWORKS
To reveal the unique properties of the system dynamics that
induced by the complex topology, we go on to investigate
the pattern formation of unsynchronizable networks by the
method of temporal phase synchronization (TPS).
A. Temporal phase synchronization
TPS is defined as follows. Let xi(t) be the time sequence
recorded on node i, we first transform it into a symbolic se-
quence θi(t) according to the following equations
θi(t) =
0, if xi(t) < 0.5,
1, if xi(t) ≥ 0.5.
Then we divide θi(t) into short segments of the equal length
n. Regarding each segment as an new element, we there-
fore have transformed the long, variable sequence xi(t) into
a short, symbolic sequence Θi(t
′). If at moment t′ we have
′) = Θj(t
′), then we say that TPS is achieved between
0.820 0.825 0.830 0.835
0 2000 4000 6000 8000 10000
10 100 1000
FIG. 3: (Color online) The on-off intermittency of the system dy-
namics nearby the LW bifurcation at ε = 0.83. (a) The time evolu-
tion of the average trajectory distance ∆X . (b) The laminar-phase
distribution of ∆X , which follows a power-law scaling with the fit-
ted exponent around 3/2. (c) The transition behavior of the average
distance 〈∆X〉 nearby the LW bifurcation point ε1, where a linear
relation is found between the two quantities.
the nodes i and j. The collection of nodes which have the
same value of Θ at moment t′ are defined as a temporarily
synchronous cluster, and all the synchronous clusters consti-
tute the temporarily pattern of the system. During the course
of system evolution, the clusters will change their shapes and
contents and the pattern will change its configuration.
In comparison with the method of complete synchroniza-
tion, the advantage we benefit from TPS is obvious: it makes
the synchronous pattern detectable. With complete synchro-
nization, it is almost impossible for two nodes to have ex-
actly the same variable at the same time. Despite the fact that
at some moments the system has already reached the high-
coherence states (formed during those contractive intervals in
Fig. 3(a)), with complete synchronization we are not able to
distinguish these states from those low-coherence ones quanti-
tatively (formed during those divergent intervals in Fig. 3(a)).
(A remedy to this difficulty seems to define the clusters by
the method of threshold truncation, i.e., nodes are regarded
as synchronized if the distance between their trajectories is
smaller than some small value. However, this definition of
synchronization will induce the problem of cluster idenfica-
tion, as the same state may generate different patterns if we
choose the different reference nodes.) On the contrary, TPS
focuses on the loose match (phase synchronization) between
the node variables over a period of time. By requiring an ex-
act match of the discrete variable Θ, the synchronous pattern
is uniquely defined; while by requiring the match of the long
sequences of θ, the ”synchronous” nodes are guaranteed with
a strong coherence.
B. Pattern evolution of the giant-cluster state
With the same set of parameters as in Fig. 3(a), by the
method of TPS we plot in Fig. 4 the time evolutions of two ba-
sic quantities of pattern evolution: the number of synchronous
clusters nc and the size of the largest cluster Lmax. It is
found that, similar to the phenomenon in complete synchro-
nization [Fig. 3(a)], on-off intermittency is also found in the
TPS quantities nc and Lmax. In Fig. 4(a) it is shown that most
of the time the system is broken into only a few number of
clusters, nc = 2 or 3, while occasionally it is broken into a
quite large number of clusters, 10 < nc < 50, or united to the
synchronization state, nc = 1. The intermittent pattern evo-
lution is also reflected on the sequence of Lmax [Fig. 4(b)],
where most of the time the size of the giant cluster is about
Lmax ≈ N , while occasionally it decreases to some small
values of Lmax < N/2. As we have discussed previously,
the main advantage we benefit from TPS is in identifying the
clusters. This advantage is clearly shown in Figs. 4(a) and
(b), where for any time instant the two quantities nc and Li
are uniquely defined. Besides cluster identification, we also
benefit from TPS in quantifying the synchronization degrees.
In specific, the different coherence states shown in Fig. 3(a)
now can be clearly quantified: high coherence states are those
of smaller nc or larger Lmax. Specially, the synchronization
state now is unambiguously defined as the moments of nc = 1
in Fig. 4(a) or, equally, the moments of Lmax = N in Fig.
4(b).
We go on to investigate the pattern evolution by statistical
analysis. The first statistic we are interested is the laminar-
phase distribution of the synchronization state, i.e. the time
intervals that nc = 1 in Fig. 4(a) or Lmax = N in Fig. 4(b).
In its original definition, laminar phase refers to the time inter-
vals τ that all node trajectories stays within a small distance
from the synchronous manifold, therefore the actual value of
τ is varying with the predefined threshold distance. This un-
certainty is overcome in TPS. As shown in Fig. 4(a), in TPS
0 1x105 2x105 3x105
0 1x105 2x105 3x105
FIG. 4: For the same set of parameters as in Fig. 3(a). The time
evolutions of the TPS quantities. (a) The number of the synchronous
clusters nc and (b) the size of the giant cluster Lmax. The synchro-
nization state is defined as the moments nc = 1 in (a) or Lmax = N
in (b).
the ”off” state refers to the moments of nc = 1 specifically.
The laminar-phase distribution of nc is plotted in Fig. 5(a). In
consistency with the distribution of complete synchronization
[Fig. 3], the laminar-phase distribution of nc also follows a
power-law scaling and has the same exponent γ ≈ −1.5±0.1.
Therefore the use of TPS, while bringing convenience to the
pattern analysis, still capture the basic properties of the on-off
intermittency. The second statistic we are interested is the size
distribution of the largest cluster, an important indicator for
the coherence degree of the system. For the Lmax sequence
plotted in Fig. 4(b), in Fig. 5(b) we plot the size distribution
of Lmax. It is seen that the probability of finding large clus-
ter Lmax ≈ N is much higher than that of small cluster of
Lmax < 500. In particular, the probability for finding clusters
of Lmax > 990 is about 20 percent and for Lmax > 990 it
is about 70 percent. Therefore, in the region of ε . ε1, the
distinct feature of the system patterns is the existence of a gi-
ant cluster. Due to this special feature, we call these states the
giant-cluster state.
Besides the giant cluster, we are also interested in the prop-
erties of the small clusters. We plot in Fig. 5(c) the distribu-
tion of nc and in Fig. 5(d) the size distribution of the small
clusters Li that surround the giant cluster in the pattern. As
shown in Fig. 5(c), the distribution of nc follows a power-law
100 101 102
0 500 1000
1 10 100
1 10 100 1000
FIG. 5: (Color online) Statistical properties of the on-off intermit-
tency plotted in Fig. 4. (a) The power-law scaling of the laminar-
phase distribution of nc. The fitted slope is about −2.3 ± 0.05. (b)
The size distribution of the size of the giant cluster. (c) The power-
law distribution of the number of small clusters nc. The fitted slope
is about −3±0.1. (d) The power-law scaling on the size distribution
of the small clusters. The fitted slope is about −1.2± 0.01.
scaling with the fixed exponent is about γ ≈ −3± 0.05. The
heterogeneous distribution of nc indicates that, in the giant-
cluster state, the system is usually broken into only a few num-
ber of clusters. An interesting finding exists in the size distri-
bution of the small clusters. As shown in Fig. 5(d), in range
Li ∈ [1, N/2] a power-law scaling is found between PLi and
Li, with the fitted exponent is about γ ≈ −1.1 ± 0.05. The
distribution of Li confirms the finding of Fig. 5(c) that the
small clusters which join or separate from the giant cluster are
usually of small size.
Combining the findings of Fig. 4 and Fig. 5, the picture
of pattern evolution in the bifurcation region ε . ε1 now be-
comes clear. Generally speaking, the evolution can be divided
into two opposite dynamical processes happening around the
giant cluster: the separation and integration of the small clus-
ters. During the separation process, the small clusters are es-
caped from the giant cluster, which weakens the dominant role
of the giant cluster and makes the pattern complicated. How-
ever, the separated clusters occupy only a small proportion of
the nodes [Fig. 5(c)], the majority nodes are still attached to
the giant cluster, which sustains the synchronization skeleton
and keeps the system on the high coherence states. At some
rare moments the giant cluster may disappears, and the pat-
tern is composed by only small clusters of Li < N/2. At
these moments, the synchronization skeleton is broken, the
pattern becomes even complicated and the system coherence
reaches its minimum. In contrast, during the process of cluster
integration, the giant cluster will increase it size by attracting
the small clusters , and gradually towards the state of global
synchronization. It should be noticed that the separation and
integration processes are uneven and are typically occurring
at the same time. For instance, during the separation process,
while the system evolution is dominated by the separation of
new small clusters from the giant cluster, there could be some
small clusters rejoin to the giant cluster.
C. Pattern evolution of the scattering-cluster state
As we further decrease the coupling strength from ε1, the
picture of pattern evolution will be totally changed. With
ε = 0.79, we plot in Fig. 6 the same statistics as in Fig. 5.
The first observation is the loss of the global synchronization
state, as can be found from the time variation of nc plotted in
Fig. 6(a). The loss of global synchronization becomes even
clear if we compare Fig. 6(a) with Fig. 4(a): in Fig. 6(a),
except the moment at t = 0, the system can never reach the
synchronization state at nc = 1 and very often it is broken into
a large number of small clusters at about nc ∼ 10
2. The fact
that the pattern is decomposed into a large number of small
clusters is also manifested by the distribution of nc, as plotted
in Fig. 6(b). Instead of the power-law distribution found in the
giant-cluster state, in the scattering-cluster state nc follows a
Gaussian distribution [Fig. 6(b)]. As ε further decreases from
ε1, the mean value of nc will shift to the larger values, as
indicated by the ε = 0.78 curve plotted in Fig. 6(b). The sec-
ond observation is the disappearance of the giant cluster. As
shown in Fig. 6(c), the size distribution of the largest cluster
also follows a Gaussian distribution, with its mean value lo-
cates at 〈Lmax〉 < N/2. The distribution of Fig. 6(c) is very
different to that of Fig. 5(b), where in Fig. 5(b) the largest (gi-
ant) cluster has size Lmax ≈ N in most of the time. As ε de-
creases, the mean value of the largest cluster 〈Lmax〉 will shift
to small values and the variance of Lmax will be decreased, as
indicated by the ε = 0.78 curve plotted in Fig. 6(c). Similar to
plot of Fig. 5(d), we have also investigated the distribution of
Li, the sizes for all the small clusters appeared in the system
evolution [Fig. 6(d)]. It is found that the distribution of Li fol-
lows a power-law distribution for Li < N/2, while having an
exponential tail for Li > N/2. Numerically we find that the
exponent of the power-law section, i.e. in rangeLi ∈ [1, 200),
is about −2± 0.05, while the fitted exponent for the exponen-
tial section is about −4.5 × 10−3 ± 2 × 10−5. These two
exponents, however, are changing with ε. As ε decreases, the
two exponents will shift to some small values.
Combining Fig. 5 and Fig. 6, we are able to outline the
transition process of network synchronization nearby the bi-
furcation points, i.e., the transition from the giant-cluster state
to the scattering-cluster state as ε leaves away from ε1. In
the region of ε . ε1, the pattern is composed by a giant
cluster and a few number of small clusters, i.e. the giant-
cluster state. As ε decreases from ε1 gradually, more and
more small clusters will be emitted out from the giant clus-
ter and, as a consequence, both the size of the giant cluster
and the fraction of synchronization time will be decreased.
Then, at about εc ≈ 0.832, the giant cluster disappears and
the pattern of the system is composed by several larger clus-
ters, of size Lmax . N/2, together with many small clusters
of heterogenous size distribution, i.e. the scattering-cluster
state. After that, as ε decreases from εc, the clusters shrink
their size by breaking into even small clusters, and the pat-
0 200 400 600 800 1000
0 100 200 300 400 500
1 10 100 1000
0.0 5.0x104 1.0x105 1.5x105 2.0x105
 g=0.78
 g=0.79
 g=0.78
 g=0.79
(d)  g=0.78
 g=0.79
FIG. 6: (Color online) The dynamical and statistical properties of
pattern evolution for ε = 0.79. (a) The time evolution of nc. (b) The
Gaussian distribution of the number of the small clusters nc. (b) The
Gaussian distribution of the size of the largest cluster Lmax. (d) The
two-segment distribution on the size of the small clusters Li. In the
region of Li < 200, Li follows a power-law distribution with fitted
exponent is about −2 ± 0.05; while for Li > 200, the distribution
is exponential with the fitted exponent is about −4.5 × 10−3 ± 2×
−5. As ε further decreases from ε1, the largest cluster becomes
even smaller and more small clusters are emitted out from it. As
illustrated by the ε = 0.78 curves plotted in (b), (c) and (d).
tern becomes even complicated. The detail transition from the
giant-cluster state to the scattering-cluster state is presented
in Fig. 7, where the average number of clusters that the sys-
tem is broken into 〈nc〉, Fig. 7(a), and the average size of
the largest cluster 〈Lmax〉, Fig. 7(b), are plotted as a func-
tion of the coupling strength in the LW bifurcation region.
The transition is found to be smooth and steady, just as we
have expected. Besides the giant cluster, another difference
between the giant-cluster and scattering-cluster states exists
in their pattern evolutions. In the giant-cluster state, while
the configuration of the giant cluster is continuously updated
by emitting or absorbing the small clusters, its main contents
are stable and do not change with time. In contrast, in the
scattering-cluster state the small clusters integrate with or sep-
arate from each other in a random fashion. Although occa-
sionally there could be some large clusters show up in the pat-
tern of the scattering-cluster state [Fig. 6(d)], these ”large”
clusters, however, are very fragile and will break into small
clusters again in a short time. This quick-dissolving prop-
erty stops the scattering-cluster state from having a high co-
herence.
VI. CHARACTERIZING THE ACTIVE NODES
In the giant-cluster state, most of the nodes are organized
into the giant cluster while few nodes, either in forms of small
group or isolated node, are separating from or joining to the
giant cluster with a high frequency. These active nodes, al-
0.79 0.80 0.81 0.82 0.83 0.84 0.85
0.79 0.80 0.81 0.82 0.83 0.84 0.85
 > (0.832,500)
FIG. 7: The transition process of the network synchronization nearby
the LW bifurcation point ε1. (a) The average number of clusters that
the system is broken into as a function of coupling strength. (b) The
average size of the largest cluster as a function of coupling strength.
Each date is an averaged result over 108 time steps.
though are few in amount, play an important role in network
synchronization. Clearly, a proper characterization of these
nodes will deepen our understandings on the system dynam-
ics and give indications to the improvement of network perfor-
mance. For instance, to improve the synchronizability of the
system, we may either remove the few most active nodes from
the network, or update their coupling strengths specifically.
In characterizing the active nodes, the following properties
are of general interest: 1) what’s the dependence of the node
activity on the network topology? can we characterize these
nodes by the known network properties such as node degree
or betweenness? 2) are their locations sensitive to the coupling
strength? and 3) what’s the effect of bifurcation type on their
locations? In the following we will explore these questions by
numerical simulations.
We first try to characterize the active nodes by their topo-
logical properties. For the giant-cluster state described in Fig.
4, we plot in Fig. 8(a) the probability pu1 that each node stays
in the giant cluster. While the majority nodes stay in the gi-
ant cluster with a high probability pu1 ≈ 1, few nodes are
of unusually small probabilities: 1 percent of the nodes have
pu1 < 0.8. One important observation of Fig. 8(a) is that
the locations of the active nodes are entangled with those of
310 320 330 340 350
0 200 400 600 800 1000
160 170 180 190 200
0 200 400 600 800 1000
 0.9514
 0.9515
 0.952(c)
 0.8342
 0.834
 0.8335
 0.83(a)
FIG. 8: (Color online) The properties of the active nodes. (a) For the
giant-cluster state shown in Fig. 4, the probability that node stays in
the giant cluster versus the node index. (b) A segment of (a) but with
different coupling strengths nearby the LW bifurcation point ε1. (c)
For the giant-cluster state (ε = 0.952) nearby the SW bifurcation
point, the probability that node stays in the giant cluster versus the
node index. (d) A segment of (c) under different coupling strengths
nearby the SW bifurcation point ε2.
the stable nodes. Noticing that in the BA growth model node
of higher index in general assume the smaller degree, the ob-
servation of Fig. 8(a) therefore indicates the independence of
the node degree to the node stability, or the inaccuracy of us-
ing degree to characterize the node activity. Specifically, in
Fig. 8(a) the 5 most unstable nodes, by a descending order of
pu1, are those of degrees k = 47, 36, 26, 10, 4, respectively.
Except the one of k = 4, all the other nodes have higher de-
grees. Another well-known topological property of complex
network is the node betweenness, which counts the number
of shortest pathes that pass through each node and actually
evaluates the node importance from the global-network point
of view. This global-network property, however, is also inca-
pable to characterize the active nodes. In Tab. 1 we list the
detail information about the 5 most active nodes in Fig. 8(a),
where the inaccuracy of node degree or node betweenness in
characterizing the active nodes are summarized.
TABLE I: For the attaching probability pi plotted in Fig. 8(a), we
list the 5 most unstable nodes and try to characterize them by a set of
topological quantities including the node index i, the attaching prob-
ability pi, the stability rank pi rank, the node degree ki, the degree
rank ki rank, the node betweenness Bi, and the betweenness rank
Bi rank.
Node index i pi pi rank ki ki rank Bi Bi rank
615 0.72797 1 5 39→537 1301 280
762 0.74424 2 5 39→537 1375 356
680 0.75416 3 4 1→338 1254 680
372 0.7591 4 6 538→645 1440 406
938 0.75972 5 4 1→338 1215 159
We go on to investigate the affection of the coupling
strength on the locations of the active nodes. In Fig. 8(b) we
fix the network topology and compare the node activities un-
der different coupling strengths nearby the bifurcation point
ε1. It is found that, despite of the changes in pu1, the lo-
cations of the active nodes are kept unchanged. That is, the
active nodes are always the first ones to escape from the gi-
ant cluster whenever the network is unsynchronizable. We
have also investigated the affection of the bifurcation type on
the locations of the active nodes. By choosing the coupling
strength nearby the SW bifurcation ε = 0.952 & ε2, we plot
in Fig. 8(c) the node attaching probability pu2 as a function
of the node index i. An interesting finding is that, comparing
to the situation of LW bifurcation [Fig. 8(a)], the locations
of the active nodes have been totally changed in Fig. (c). In
Tab. 2 we list the detail information about the 5 most active
nodes in Fig. 8(c), again their locations can not be predicted
by the node degree or betweenness. Similar to the LW bifur-
cation, the locations of the active nodes are also independent
to the coupling strength at the SW bifurcation, as shown in
Fig. 8(d).
TABLE II: Similar to Tab. I but for the attaching probability pi plot-
ted in Fig. 8(c). Comparing to Tab. I, one important observation is
the changed locations of the active nodes due to the changed bifurca-
tion type.
Node index i pi pi rank ki ki rank Bi Bi rank
43 0.78196 1 9 779→813 2847 813
35 0.78969 2 18 936→940 8513 953
714 0.795 3 4 1→338 1215 158
130 0.79652 4 13 886→901 4200 886
154 0.19944 5 10 814→846 2898 815
Previous studies about network synchronization have
shown that, while individually it is difficult to predict the dy-
namical behavior of each node, the average performance of
an ensemble of nodes of the same network properties do have
some reliable characters. For instance, it has been shown that
in complex networks the high-degree nodes are on average
more synchronizable than the low-degree ones [19]. Regard-
ing to the problem of node activities, it is natural to ask the
similar question: are the high-degree nodes more synchro-
nized than the low-degree nodes? In Fig. 9 we plot the av-
erage attaching probability 〈pu1〉k as a function of degree k.
Still, we can not find a clear dependence of 〈pu1〉k on k.
VII. DISCUSSIONS AND CONCLUSION
It is worthy of note that our studies of active nodes are only
focused on the giant-cluster state, and the purpose is to under-
stand their dynamics and reveal their properties. By ensemble
0 30 60 90 120
FIG. 9: (Color online) The average attaching probability 〈pu1〉k as
a function of node degree k. On average, the 5 most unstable nodes
are those of degrees k = 47, 36, 26, 10, 4. Still, we can not find a
clear dependence between 〈pu1〉k and k.
average, we may able to improve our prediction of the active
nodes, say for example the dependence of 〈pu1〉k on k in Fig.
9 may be smoothed if we average the results over a large num-
ber of network realizations. Such an improvement, however,
comes at the cost of the decreased prediction accuracy due to
the increased candidates. Taking Fig. 9 as an example, al-
though it is noticed that nodes of k = 4 in general are more
active than those of other degrees, only one of them is listed
as the 5 most unstable nodes [Tab. 1]. In specific, among the
total number of 338 nodes which have degree k = 4, most
of them are tightly attracted to the giant cluster (90 percent of
them have attaching probabilities pu1 > 0.95). Therefore, in
terms of precise predication, the average method is infeasible
in practice.
Beside node degree and betweenness, we have also checked
the dependence of the property of node activity to some other
well-known network properties such as the clustering coeffi-
cient, the modularity, and the assortativity. However, none of
them is suitable to characterize the active nodes, their perfor-
mance is very similar to that of the node degree described in
Tab. 1 and Tab. 2. Our study thus suggests that, to give a pre-
cise prediction to the active nodes, we may need to develop
some new quantities.
Despite of the amount of studies carried on network syn-
chronization, to the best of our knowledge, we are the first to
study the nonstationary pattern in unsynchronizable complex
networks. In Ref. [9] the authors have discussed the transient
process of global synchronization in complex networks, but
their study are concentrating on the synchronizable state in
which, during the course of system evolution, small clusters
integrate into larger clusters monotonically and finally reach
the synchronization state. After that, the system will always
stay on the synchronization state. Our works are also different
to the studies of Refs. [10, 20]. Similar to our works, in these
studies the authors also consider the problem of pattern for-
mation in unsynchronizable networks, but their interests are
focused on the stationary pattern of the system. That is, the
size and contents of the clusters do not change with time. In
contrast, in our studies both the size and contents of the clus-
ters are variable.
In summary, we have reported and investigated a kind of
new phenomena in network synchronization: the nonstation-
ary pattern. That is, the final state of the network settles nei-
ther to the synchronization state nor to any stationary state
of fixed pattern, the system is travelling among all the possi-
ble patterns in an intermittent fashion (the pattern can be of
any configuration, but its probability of showing up is pattern-
dependent). We attribute this nonstationarity to the asymmet-
ric topology of the complex networks, and its dynamical ori-
gin can be understood from the property of the finite-time
Lyapunov exponent associated to the desynchronized mode.
Two types of synchronization formats, the complete synchro-
nization and the temporal phase synchronization, have been
employed to detect the nonstationary dynamics. For coupling
strength immediately out of the stable region, the pattern evo-
lution is characterized by the process of on-off intermittency
and the existence of the giant-cluster; while if the coupling
strength is far away from the bifurcation points, the pattern
evolution is signatured by the random interactivities among
the number of small clusters. A remarkable finding is that,
in the giant-cluster state the locations of the active nodes are
independent of the coupling strength but are sensitive to the
bifurcation types. The active nodes, however, can not be char-
acterized by the currently known network properties, further
investigations about their identification are necessary. While
we are hoping our studies about nonstationary pattern could
give some new understandings to the dynamics of coupled
complex systems, we also hope that our findings about unsyn-
chronizable networks could be used to some practical prob-
lems where system maintains their normal functions only un-
der the unsynchronizable states, for example the problem of
epileptic seizers [21].
[1] D.J. Watts and S.H. Strogatz, Nature 393, 440 (1998); A.-L.
Barabási and R. Albert, Science 286, 509 (1999); R. Albert and
A.-L. Barabási, Rev. Mod. Phys. 74, 47 (2002).
[2] S. Boccaletti and L.M. Pecora, Chaos 16, 015101 (2006); A.
E. Motter, M. A. Matı́as, J. Kurths, E. Ott, Physica D 224, 7
(2006); S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and
D.-U. Hwang, Phys. Rep. 424, 175 (2006).
[3] X.F. Wang and G. Chen, Int. J. Bifurcation Chaos Appl. Sci.
Eng. 12, 187 (2002).
[4] T. Nishikawa, A. E. Motter, Y.-C. Lai, and F. C. Hoppensteadt,
Phys. Rev. Lett. 91, 014101 (2003).
[5] A.E. Motter, C. Zhou, and J. Kurths, Europohys. Lett. 69, 334
(2005); Phys. Rev. E 71, 016116 (2005); AIP Conf. Proc. 776,
201 (2005); C. Zhou, A.E. Motter, and J. Kurths, Phys. Rev.
Lett. 96, 034101 (2006).
[6] M. Chavez, D.-U. Hwang, A. Amann, H.G.E. Hentschel, and S.
Boccaletti, Phys. Rev. Lett. 94, 218701 (2005); D.-U. Hwang,
M. Chavez, A.Amann, and S. Boccaletti, Phys. Rev. Lett. 94,
138701 (2005).
[7] L.M. Pecora and T.L. Carroll, Phys. Rev. Lett. 80, 2109 (1998);
M. Barahona and L. M. Pecora, ibid, 89, 054101 (2002).
[8] D.-S. Lee, Phys. Rev. E 72, 026208 (2005); J. Gómez-
Gardeñes, Y. Moreno, and A. Arenas, Phys. Rev. Lett. 98,
034101 (2007).
[9] A. Arenas, A. Dı́az-Guilera, and C. J. Pérez-Vicente, Phys. Rev.
Lett. 96, 114102 (2006); C. Zhou, L. Zemanová, G. Zamora, C.
C. Hilgetag, and J. Kurths, ibid, 97, 238103 (2006)
[10] S. Boccaletti, M. Ivancheko, V. Latora, A. Pluchino, and A.
Rapisarda, Preprint physics/0607179 (2006).
[11] M. Zhan, Z.G. Zheng, G. Hu, and X.H. Peng, Phys. Rev. E
62, 3552 (2000); Y. Zhang, G. Hu, H.A. Cerdeira, S. Chen, T.
Braun, and Y. Yao, Phys. Rev. E 63, 026211 (2001); B. Ao and
Z. Zheng, Europhys. Lett. 74, 229 (2006); X. Zhang, M. Fu, J.
Xiao, and G. Hu, Phys. Rev. E 74, 015202 (2006).
[12] E. Ott and J.C. Sommerer, Phys. Lett. A 188, 39 (1994); M.
Ding and W. Yang, Phys. Rev. E 56, 4009 (1997); S.H. Wang,
J. Xiao, X.G. Wang, B. Hu, and G. Hu, Eur. Phys. J. B 30, 571
(2002); J.G. Restrepo, E. Ott, and B. Hunt, Phys. Rev. Lett. 93,
114101 (2004).
[13] X.G. Wang, Y.-C. Lai, and C.-H. Lai, Preprint nlin.CD/0608035
(2006).
[14] L. A. Bunimovich, A. Lambert, and R. Lima, J. Stat. Phys. 65,
253 (1990); J.F. Heagy, T.L. Carroll, and L.M. Pecora, Phys.
Rev. Lett. 73, 3528 (1994); ibid, 74, 4185 (1995).
[15] M.A. Matias, V.P. Munuzuri, M.N. Lorenzo, I.P. Marino, and
V.P. Villar, Phys. Rev. Lett. 78, 219 (1997); G. Hu, J. Yang, and
W. Liu, Phys. Rev. E 58, 4440 (1998).
[16] A. Pikovsky and U. Feudel, Chaos 5, 253 (1995); X. Wang, M.
Zhan, C. H. Lai, and Y.-C. Lai, Phys. Rev. Lett. 92, 074102
(2004).
[17] N. Platt, E.A. Spiegel, and C. Tresser, Phys. Rev. Lett. 70,
279 (1993); E. Ott and J. C. Sommerer, Phys. Lett. A 188, 39
(1994); Y. Nagai and Y.-C. Lai, Phys. Rev. E 55, 1251 (1997).
[18] A.S. Pikovsky, M.G. Rosenblum, and J. Kurths, Synchroniza-
tion: A Universal Concept in Nonlinear Science (Cambridge
University Press, Cambridge, UK, 2001); S. Strogatz, Sync:
The Emerging Science of Spontaneous Order (Hyperion, New
York, 2003).
[19] C. Zhou and J. Kurths, Chaos 16, 015104 (2006).
[20] S. Jalan, R.E. Amritkar, and C.-K. Hu, Phys. Rev. E 72, 016211
(2005); ibid, 72, 016212 (2005).
[21] Y.-C. Lai, M. G. Frei, I. Osorio, and L. Huang, Phys. Rev. Letts.
92, 108102 (2007).
http://arxiv.org/abs/physics/0607179
http://arxiv.org/abs/nlin/0608035
