Many realistic epidemic networks display statistically synchronous behavior which we will refer to as epidemic synchronization. However, to the best of our knowledge, there has been no theoretical study of epidemic synchronization. In fact, in many cases, synchronization and epidemic behavior can arise simultaneously and interplay adaptively. In this paper, we first construct mathematical models of epidemic synchronization, based on traditional dynamical models on complex networks, by applying the adaptive mechanisms observed in real networks. Then, we study the relationship between the epidemic rate and synchronization stability of these models and, in particular, obtain the conditions of local and global stability for epidemic synchronization. Finally, we perform numerical analysis to verify our theoretical results. This work is the first to draw a theoretical bridge between epidemic transmission and synchronization dynamics and will be beneficial to the study of control and the analysis of the epidemics on complex networks.

Exploring the correlation between different dynamical behaviors which may appear in complex networks is an important and interesting pursuit on the way toward understanding them better. Among these different dynamical behaviors, synchronization and epidemic spreading on networks are closely related. For example, with the spread of an infective disease, people may reduce the frequency of inter-personal contact and take collective protective measures with increased regularity (for example, by washing hands frequently with clear water and soap, avoiding going to crowded place, and so on). This means that the spread of information between people during transmission of an epidemic disease can induce spontaneously collective risk-minimisation behaviour in spite of (or at least independently of) the actual disease pathology. Similar phenomena can also be found in many animal-borne infections. Hence, synchronization of individual’s behavior and their epidemic behavior on networks can occur simultaneously and interplay adaptively. However, very little theoretical work has been done to consider these two dynamical behaviors together. In this work, we address this deficit. We will investigate mathematically the correlation between the dynamical synchronization and the epidemic behavior on complex networks for us to understand them in depth.

It is well known that complex networks can accurately describe the topological structure of many dynamical systems in the real world. For example, the world wide web can be considered as a network of websites with hyperlinks. The Internet is formed by routers with physical connections. Similarly, the human brain can be also regarded as a system of neurons linked by synapses. Many other examples can be also found: the global economy, market, food cycle, metabolism, disease spreading, computer virus, energy diffusion, etc. Recently, research in the area of complex network has advanced significantly and will lead to a deeper study and more practical applications in future. Synchronization dynamics and epidemic dynamics are two of the main research fields in complex network science. In the last ten years, with the important discovery of small-world and scale-free networks, synchronization in complex networks has attracted great attention.1–13 Epidemic transmission in complex networks has also become an important research focus as it offers potential to advance the simple assumptions of homogeneously mixed models.14–19 

Since synchronization and epidemic spreading are apparently two quite distinct behaviors, there is very little work to consider them together. Actually, synchronization and epidemic behavior on networks can arise simultaneously and interplay adaptively. Many realistic epidemic networks have displayed statistically synchronized behavior that we define as the epidemic synchronization which, for the moment, can be classified as epidemic induced synchronization and epidemic characteristic synchronization. For examples, with the spread of SARS (the severe acute respiratory syndrome), Bird Flu, or H1N1 influenza, people will spontaneously take some protective measures such as washing hands frequently with clear water and soap, avoiding going to crowded place, and so on in order to improve self-protection.20 When rhinoceroses and hippopotami are infected with epidemic ringworm, they will spontaneously go to a lakefront to have a bath in mud which will also provide either cure or protection. As the spread of these epidemic diseases becomes weak, these spontaneous protective behaviors (in either humans or animals) will decrease accordingly. Here, we choose to describe these consistent human or animal behaviors, which result from epidemic spreading, as the epidemic induced synchronization. For another example, weekly measles case reports21 for Birmingham, Newcastle, Cambridge, and Norwich between 1944 and 1958 have shown their singular oscillatory behavior and an in-phase synchronized pattern between Birmingham and Newcastle, but an anti-phase pattern between Cambridge and Norwich. Similarly, the reported cases of syphilis and gonorrhoea for the Midwest, South, Northeast, and West cities in United States between 1941 and 2002 have also exhibited the oscillatory characteristic and phase synchronized patterns.22 In addition, collective synchronization induced by epidemic dynamics on complex networks has also been studied numerically.23 For integrality, we define here these consistent oscillation of measures, which quantify the epidemic degree or other characteristics and do not refer to concrete human or animal behaviors, in epidemic networks as the epidemic characteristic synchronization. In this paper, we just address the epidemic induced synchronization.

The spatial synchronization of epidemics is considered to be a key factor affecting their long-term patterns of persistence and extinction.24 However, there are very few theoretical works exploring the real reason for epidemic synchronization apart from some numerical analyses from the view of statistics. So, it is significant to seek a reliable theory for epidemic synchronization. To this end, in this paper, we first establish the mathematical models to realize the epidemic synchronization and make further analysis of these models to reveal their exhibition mechanisms, based on the theory of dynamical systems.

This paper is organized as follows. In Sec. II, we introduce the traditional models of complex dynamical network and the epidemic network. In Sec. III, the models of SIS (the susceptible-infected-susceptible) epidemic synchronization and SIR (the susceptible-infected-recovered or the susceptible-infected-removed) epidemic synchronization are constructed. Then, we investigate the local and global conditions of epidemic synchronization with respect to the epidemic rate of these models. Finally, in Sec. IV, we use some numerical simulations to verify our theoretical results obtained in Sec. III. Section V concludes this paper and lists some possible further works.

Without loss of generality, we introduce a complex dynamical network5 with linear coupling which can be described as

x·i(t)=f(xi(t))+cj=1NaijHxj(t),i=1,2,,N,
(1)

where xi(t)Rn denotes the state variable of the ith node at the time t, and the function f(·) defines the local dynamics of each node and is supposed to be chaotic. The constant c > 0 is coupling strength and the matrix HRn×n represents the inner-coupling matrix which is a constant 0 – 1 matrix linking coupled variables, and we assume it is positive. The coupling matrix A = (aij)N × N with zero-sum rows shows the coupling configuration of the network. If nodes i and j are connected, then aij = aji = 1, otherwise aij = aji = 0. The diagonal elements of the coupling matrix A are

aii=-j=1,jiNaij=-ki,i=1,2,,N,
(2)

where ki denotes the degree of node i. With these assumptions, the eigenvalues of matrix A can be given by 0=λ1λ2λN. Hence, by matrix theory, there exists unitary matrix U such that A = UΛUT, where UTU = I, Λ = diag(λ1, λ2, … ,λN).

Now, the standard SIS and SIR epidemic models17,19 in a complex network can be written as

I·k(t)=λk[1-Ik(t)]Θ(t,k)-Ik(t),k=1,2,,dm
(3)

and

{S·k(t)=-λkSk(t)Θ(t,k),I·k(t)=λkSk(t)Θ(t,k)-Ik(t),R·k(t)=Ik(t),k=1,2,,dm,
(4)

respectively. Here, the epidemic rate λ(0,1] denotes the probability with which each susceptible node is infected if it is connected to one infected node and dm represents the maximal degree in the network. The variables Sk(t), Ik(t), and Rk(t) denote the densities of susceptible, infected, and removed nodes (individuals) with connectivity (contact) k at time t, respectively. These variables satisfy the normalization condition Sk(t) + Ik(t) = 1 or Sk(t) + Ik(t) + Rk(t) = 1 for all k-classes in model (3) or (4). The term Θ(t,k) gives the probability that a randomly chosen link emanating from a node of connectivity k leads to an infected node. Moreover, Θ(t,k) has the form

Θ(t,k)=k'p(k'|k)Ik'(t),
(5)

where the conditional probability p(k'|k) means that a randomly chosen link emanating from a node of connectivity k leads to a node of connectivity k′. We suppose that the connectivities of nodes in the whole network are uncorrelated, i.e., p(k'|k)=k'p(k')/k, where k=ssp(s). The further details of models (3) and (4) can be found in Refs. 17 and 19.

Based on the above traditional models (1) and (3), we can construct the following model of SIS epidemic synchronization:

{ẋi(t)=f(xi(t))+c(t)j=1NaijHxj(t),İk(t)=λk[1Ik(t)]Θ(t,k)Ik(t),ċ(t)=αI(t)Nj=1Ns(t)xj(t)21+s(t)xj(t)2,
(6)

where i = 1, 2, …, N, k = 1, 2, … , dm, parameter α > 0, and I(t)=k=1dmp(k)Ik(t). The initial condition of system (6) can be set as follows. The initial state xi(0) is chosen randomly from real number set and Ik(0) = ρ, c(0) = 0 with 0<ρ1. The state variable s(t) is the synchronous state of system (6) with respect to epidemic induced synchronization or epidemic characteristic synchronization.

By the similar entrainment mechanism between the above traditional models (1) and (4), we can design the following model of SIR epidemic synchronization as

{ẋi(t)=f(xi(t))+c(t)j=1NaijHxj(t),Ṡk(t)=λkSk(t)Θ(t,k),İk(t)=λkSk(t)Θ(t,k)Ik(t),Ṙk(t)=Ik(t),ċ(t)=αI(t)Nj=1Ns(t)xj(t)21+s(t)xj(t)2,
(7)

where i = 1, 2, … , N, k = 1, 2, … , dm, parameter α > 0, and epidemic prevalence I(t)=k=1dmp(k)Ik(t).The initial condition of system (7) can be set as follows. The initial state xi(0) is chosen randomly from real number set and Sk(0) = ρ, Ik(0) = 1 – ρ, Rk(0) = 0, c(0) = 0, where 0<ρ1.

When we consider the epidemic induced synchronization, xi(t) indicates the state variable of node-i in the epidemic network. In this case, the models (1) and (3) (similarly for Eqs. (1) and (4)) have the same topological structure, i.e., dynamical process of individuals and epidemic process arise simultaneously in the same network. And f defines the locally dynamical behavior which indicates the summation of various behavior of individual node concerning epidemic. By supposing that f is chaotic, we mean that the functional behaviour of each individual node is active, dynamic, and sufficiently complex.

When we take the epidemic characteristic synchronization into account, the variables xi(t) denotes the various factors influencing the epidemic transmission (such as the weekly infection rate of cities) and f shows the locally dynamical behavior of that variable. In this case, the models (1) and (3) (similarly for Eqs. (1) and (4)) may have different general topological structures. In this paper, we just address the first case where the synchronization is induced by the epidemic dynamics, which is realized by the adaptive coupling strength c(t) mainly depending on the change of I(t) and synchronization error.

In the third equation of model (6), the change of c(t) is controlled synthetically by I(t) and j=1Ns(t)-xj(t)21+s(t)-xj(t)2, where I(t) denotes the average density of infected nodes and this summation term can measure synchronization degree in a network. On one hand, individuals will send spontaneously safeguard information more frequently to protect themselves when disease prevalence becomes larger. So, it is reasonable to conclude that the rate of change of the coupling strength c·(t) is directly proportional to the density I(t) of infected nodes. On the other hand, when the collective protective behavior increases sufficiently, communication of safeguard information among individuals will become stable since they have come to an agreement of protection. Thus, the proportional relation between c·(t) and synchronization error j=1Ns(t)-xj(t)21+s(t)-xj(t)2 is also valid. We take the same consideration in the fourth equation of model (7).

For the network system (1), we have the following basic result.

Lemma 1: Considering the network (1) with chaotic individual nodes, the coupling strength c = c(t), H = IN, and the maximal Lyapunov Exponent hmax of function f, if there is T > 0 such that

c(t)>hmax|λ2|,fort>T,
(8)

then the synchronization of network (1) is exponentially stable.

Proof: By reference,25 we know that the n(N – 1) Lyapunov exponents of network (1) on the synchronization manifold can be expressed as μi(λk)=hi+limsuptc(t)λk, where hi is the Lyapunov exponent of f and i = 1, 2, … , n, k = 2, 3, … , N. If there is T > 0 such that c(t)>hmax|λ2| for all t > T, then we get hi+limsuptc(t)λk<0 for all i = 1, 2, … , n, k = 2, 3, … , N. This means all n(N – 1) Lyapunov exponents of network (1) on the synchronization manifold are negative. According to the transverse stability of synchronization,25 we obtain the exponential stability of network (1).

Theorem 1: Suppose that λc > 0 is the epidemic threshold of system (3) and f(·) is a chaotic function of the corresponding model (6) with H = IN. If the epidemic rate λ > λc in system (6), then xi(t), i = 1, 2, … ,N can achieve synchronization.

Proof: If λ > λc, by the definition of epidemic threshold, there is a constant I*(0,1] such that limt+I(t)=I*. By defining δ(t)=αI(t)Nj=1Ns(t)-xj(t)21+s(t)-xj(t)2, we know δ(t) ≥ 0. Now, we show that xi(t), i = 1, 2, … ,N can realize synchronization in model (6).

If this is not the case, we obviously have limt+δ(t)0. Now, we will prove limtc(t)=+. In order to verify this limitation, we first show c(t) is boundless. Otherwise, there exists M > 0 such that |c(t)| < M for all t ≥ 0. By noting c·(t)=δ(t)0, we can deduce that c(t) must have limitation with t → + , which means limtδ(t)=0. So c(t) is unbounded in (0, + ). By combining the monotone property of c(t), we get limtc(t)=+. On the other hand, by Lemma 1, xi(t), i = 1, 2, … , N can achieve synchronization if c(t)>hmax/|λ2|. This contradiction concludes this theorem.

Theorem 2: Suppose that λc > 0 is the epidemic threshold of system (3) and f(·) is a chaotic function of the corresponding model (7) with H = IN. If the epidemic rate λ > λc in system (7), then xi(t), i = 1, 2, … , N can achieve synchronization.

Proof: The proof is similar to the proof of theorem 1, we omit it here for simplicity.

Now, we consider the infinite state of coupling strength c(t).

Theorem 3: Considering the model (6) or (7), for every λ(0,1), there is always a constant c*>0, such that limtc(t)=c*.

Proof: Obviously, for the model of epidemic synchronization (6) or (7), its epidemic threshold satisfies λc(0,1) in its corresponding epidemic model. Now, we will prove this theorem for two cases with respect to λ.

On one hand, if 1 < λλc, by the above analysis, we know that xi(t), i = 1, 2, … ,N can realize synchronization. This means that limts(t)-xj(t)=0 for every j{1,2,,N}. Moreover, it is easy to see that 0<αI(t)NαN. So, we get limtc·(t)=0 for the model of SIS epidemic synchronization (6) or SIR epidemic synchronization (7). By combining the initial condition c(0)=0,c·(0)>0 and inequality c·(t)0, we find that there is always a constant c* > 0, such that limtc(t)=c*.

On the other hand, if 0 < λλc, we have limtI(t)=0 for the models (6) and (7) by the definition of epidemic threshold. In addition, it is obvious to see that

αNj=1Ns(t)xj(t)21+s(t)xj(t)2<α.

Thus, for the model of SIS epidemic synchronization (6) or SIR epidemic synchronization (7), we obtain limtc·(t)=0. In this case, by combining the initial condition c(0)=0, c·(0)>0, and inequality c·(t)0, we obtain the same result. Hence, the Theorem is proved.

Based on the paper,10 the synchronous state in network (6) or (7) can be defined as s(t)=1Ni=1Nxi(t). Hence, the corresponding synchronization errors can be set as ei(t) = xi(t) – s(t), i = 1, 2, … ,N. It is easy to obtain i=1Nei(t)=0 and s·(t)=1Ni=1Nf(xi(t)). Consequently, the error system can be described as

e·i(t)=f(xi(t))-f(s(t))+c(t)j=1NaijHej(t)+g(t),i=1,2,,N,
(9)

where g(t)=f(s(t))-1Ni=1Nf(xi(t)).

Hypothesis 1: Suppose that P = diag(p1,p2, … ,pn) is a positive matrix. If there is a constant ξ, such that for all x(t),y(t)Rn, and t > 0, then we always have that

(x-y)TP[f(x)-f(y)]ξ(x-y)T(x-y).
(10)

By letting F(t)=(f(x1(t))Tf(s(t))T,,f(x1(t))Tf(s(t))T)T, G(t)=(gT,,gT)T, and e(t)=(e1T,,eNT)T, then the system (9) is rewritten as

e·(t)=F(t)+c(t)(AH)e(t)+(INIn)G(t),
(11)

where is Kronecker product.

Theorem 4: Suppose that λc > 0 is the epidemic threshold of system (3). If λ > λc in system (6), then the synchronous manifold of system (6) is globally asymptotically stable.

Proof: We construct the following Lyapunov function candidate

V(t)=12eT(t)(INP)e(t)+12β(c0-c(t))2,

where β  = − λ2λmin(PH) > 0, λmin(PH) denote the minimal eigenvalue of matrix PH and c0 is a undetermined constant. The derivative of V(t) with respect to t along the solution of system (6) is given by

dV(t)dt=i=1Nei(t)TP[f(xi(t))-f(s(t))]+c(t)e(t)T(APH)e(t)+e(t)T(INP)G(t)-β(c0-c(t))c·(t)=i=1Nei(t)TP[f(xi(t))-f(s(t))]+c(t)e(t)T(APH)e(t)-β(c0-c(t))c·(t)ξi=1Nei(t)Tei(t)+c(t)e(t)T(APH)e(t)-β(c0-c(t))c·(t).
(12)

Now, introducing a transformation y(t)=(y1T(t),,yNT(t))T=(UTIn)e(t) and combining Eq. (12), we obtain

dV(t)dtξi=1Nei(t)Tei(t)+c(t)yT(t)(ΛPH)y(t)-β(c0-c(t))c·(t)=ξi=1Nei(t)Tei(t)-βc0c·(t)+c(t)yT(t)(ΛPH)y(t)+βc(t)c·(t)ξi=1Nei(t)Tei(t)+c0λ2λmin(PH)αI(t)Ni=1Nei(t)Tei(t)1+ei(t)Tei(t)+c(t)λ2λmin(PH)αI(t)Ni=1Nei(t)Tei(t)1+ei(t)Tei(t)-c(t)λ2λmin(PH)c·(t).
(13)

If epidemic rate λ > λc, then there exists I*(0,1], such that limt+I(t)=I*. By choosing ɛ0(0,I*), then there is t0 such that I(t)>I*-ɛ0>0 for all t > t0. When t > t0, by Eq. (13) we have

dV(t)dtξi=1Nei(t)Tei(t)+c0λ2λmin(PH)α(I*-ɛ0)Ni=1Nei(t)Tei(t)1+ei(t)Tei(t)+c(t)λ2λmin(PH)[αI(t)Ni=1Nei(t)Tei(t)1+ei(t)Tei(t)-c·(t)]=ξi=1Nei(t)Tei(t)+c0αλ2λmin(PH)(I*-ɛ0)Ni=1Nei(t)Tei(t)1+ei(t)Tei(t).

Thus, we can select an adequately large constant c0, such that dV(t)dt0. On the other hand, it is easy to see that λmin(P)2eT(t)e(t)+12β(c0-c(t))2V(t)λmax(P)2eT(t)e(t)+12β(c0-c(t))2. So, the Lyapunov function V(t) has infinitesimal upper bound and is infinitely large.26,27 Therefore, the zero solution of error system (9) is globally asymptotically stable, i.e., the synchronisation manifold of system (6) is globally asymptotically stable.

Similarly, if epidemic rate λ > λc, the synchronisation manifold of system (7) is also globally asymptotically stable.

Remark 1: Actually, we just need the existence of the above constant ξ in Hypothesis 1 to achieve the global stability of epidemic synchronization, as shown by the above analysis. So, it is convenient to use Theorem 4 for analysis in realistic applications.

Remark 2:. In fact, since some certain synchronization patterns can accelerate or weaken epidemic prevalence, it is more realistic to consider the effect of these synchronization dynamics on epidemics synchronously, as we investigate the relationship between the epidemic dynamics and synchronization dynamics on a network. If we embed this anti-effect function in model (6) or (7), it must influence the epidemic process, i.e., either enhancing or reducing the epidemic prevalence I(t). But from the results obtained above, our synchronization conditions are all independent on this epidemic process, but depend only on the epidemic rate. So, the above synchronization conditions are still valid in this case.

To verify the above results, we now investigate numerically the model of SIS epidemic synchronization (6). Since this analysis can be generalized analogously to the model of SIR epidemic synchronization (7), we omit this generalization for simplicity. The network embedded in model (6) is made to be the BA (the Barabasi and Albert model) scale-free (preferential attachment) network28 with size N = 500. This network evolves from initial network with size m0 = 4 and we add each new node with m = 3 new edges. Without loss of generality, we suppose that the f in model (6) is defined as the chaotic Lorenz oscillation. Certainly, this assumption is rather difficult to justify from the point of view of a physical epidemic transmission process, but we just use it for numerical simulations. This oscillation can be described as

{dx1(t)dt=a1(x2(t)-x1(t))dx2(t)dt=a2x1(t)-x2(t)-x1(t)x3(t)dx3(t)dt=x1(t)x2(t)-a3x3(t)
(14)

where parameters a1 = 10, a2 = 28, and a3 = (8/3). H is chosen as the identity matrix. The Figures 1, 2, and 3 show the changes in synchronization error e(t), epidemic prevalence I(t), and coupling strength c(t) in the model (6) with epidemic rate λ = 0.02, 0.1, and 0.2, respectively. Here, the synchronization error is defined as e(t)=j=2Nx1(t)-xj(t)/N. From Fig 1, we can see that I(t) → 0(t), i.e., the epidemic does not break out. In this case, the synchronization error e(t) does not converge to zero, which means that the epidemic dynamics cannot successfully induce the synchronization of individuals under small epidemic rate λ = 0.02. From Fig 2, we can see that I(t) converges to a positive number and e(t) converges to zero, which implies that the epidemic dynamics can induce successfully the synchronization of individuals with larger epidemic rate λ = 0.1. Increasing further the epidemic rate to 0.2, we find that the epidemic dynamics can not only induce successfully the synchronization, but also enhance the speed of synchronization. These simulations show that if the epidemic transmission spreads more easily on the network, then the epidemic dynamics can induce the synchronization of individuals more effectively. This is consistent with the attained theoretic results in Sec. IIIan aggressive disease leads to a strong response from individuals.

FIG. 1.

The changes of synchronization error e(t), epidemic prevalence I(t), and coupling strength c(t) in model (6) with epidemic rate λ = 0.02. The network is BA scale-free network with size N = 500, the maximal degree dm = 60 and 192 nodes with the minimal degree 3. Other parameters are set as α = N – 1, initial infection density ρ3 = 0.03125, ρk = 0, k = 4,…,dm, i.e., there are originally 6 infected nodes with degree 3 and the other nodes are all susceptible.

FIG. 1.

The changes of synchronization error e(t), epidemic prevalence I(t), and coupling strength c(t) in model (6) with epidemic rate λ = 0.02. The network is BA scale-free network with size N = 500, the maximal degree dm = 60 and 192 nodes with the minimal degree 3. Other parameters are set as α = N – 1, initial infection density ρ3 = 0.03125, ρk = 0, k = 4,…,dm, i.e., there are originally 6 infected nodes with degree 3 and the other nodes are all susceptible.

Close modal
FIG. 2.

The changes of synchronization error e(t), epidemic prevalence I(t), and coupling strength c(t) in model (6) with epidemic rate λ = 0.1. The network is BA scale-free network with size N = 500, the maximal degree dm = 60 and 192 nodes with the minimal degree 3. Other parameters are set as α = N – 1, initial infection density ρ3 = 0.03125, ρk = 0, k = 4,…,dm, i.e., there are originally 6 infected nodes with degree 3 and the other nodes are all susceptible.

FIG. 2.

The changes of synchronization error e(t), epidemic prevalence I(t), and coupling strength c(t) in model (6) with epidemic rate λ = 0.1. The network is BA scale-free network with size N = 500, the maximal degree dm = 60 and 192 nodes with the minimal degree 3. Other parameters are set as α = N – 1, initial infection density ρ3 = 0.03125, ρk = 0, k = 4,…,dm, i.e., there are originally 6 infected nodes with degree 3 and the other nodes are all susceptible.

Close modal
FIG. 3.

The changes of synchronization error e(t), epidemic prevalence I(t), and coupling strength c(t) in model (6) with epidemic rate λ = 0.2. The network is BA scale-free network with size N = 500, the maximal degree dm = 60 and 192 nodes with the minimal degree 3. Other parameters are set as α = N – 1, initial infection density ρ3 = 0.03125, ρk = 0, k = 4,…,dm, i.e., there are originally 6 infected nodes with degree 3 and the other nodes are all susceptible.

FIG. 3.

The changes of synchronization error e(t), epidemic prevalence I(t), and coupling strength c(t) in model (6) with epidemic rate λ = 0.2. The network is BA scale-free network with size N = 500, the maximal degree dm = 60 and 192 nodes with the minimal degree 3. Other parameters are set as α = N – 1, initial infection density ρ3 = 0.03125, ρk = 0, k = 4,…,dm, i.e., there are originally 6 infected nodes with degree 3 and the other nodes are all susceptible.

Close modal

Based on the synchronization behavior of individuals induced by epidemic dynamics in real epidemic networks, this paper has proposed some mathematical models of epidemic synchronization which can characterize this kind of phenomenon very well. We have studied the relationship between the epidemic rate and synchronization stability of these models and obtained a very explicit condition for synchronization with respect to the epidemic rate, i.e., λ > λc. Numerical simulations show that if the epidemic disease breaks out more easily, then the epidemic dynamics can induce the synchronization of its individuals more effectively. This conclusion accords very well with the characteristics of real epidemic networks.

In future, we will extend our work to address the following two issues. First, since in realistic epidemic networks the individuals are generally different, the local dynamics of nodes can not be identical in a whole network. So, it must be more appropriate to consider network models with both community structure and (distinct) individual behavior. Second, we anticipate that some certain synchronization patterns can accelerate or weaken epidemic prevalence. So, it may be very significant to consider the effect of these synchronization dynamics on epidemics when we study the epidemic synchronization on complex networks.

This research was supported jointly by National Natural Science Foundation of China (Nos. 61004101, 60964006, 11072136, and 11061010). The authors gratefully acknowledge the support of the Guangxi Natural Science Foundation Program (Nos. 2011GXNSFB018059, 2011GXNSFA018136), Shanghai Leading Academic Discipline Project (No. S30104), and a grant from Guilin University of Electronic Technology.

1.
C. W.
Wu
and
L. O.
Chua
,
IEEE Trans. Circ. Syst., I: Fundam. Theory Appl.
42
,
430
(
1995
).
2.
L. M.
Pecora
and
T. L.
Carroll
,
Phys. Rev. Lett.
80
,
2109
(
1998
).
3.
I. V.
Belykh
,
V. N.
Belykh
, and
M.
Hasler
,
Physica D
195
,
159
(
2004
).
4.
A.
Pikovsky
,
M.
Rosenblum
, and
J.
Kurths
,
Synchronization: A Universal Concept in Nonlinear Science
(
Cambridge University Press
,
Cambridge
,
2001
).
5.
X. F.
Wang
and
G. R.
Chen
,
IEEE Trans. Circuits Syst., I: Fundam. Theory Appl.
49
,
54
(
2002
).
6.
Z. J.
Ma
,
Z. R.
Liu
, and
G.
Zhang
,
Chaos
16
,
023103
(
2006
).
7.
W. W.
Yu
,
J. D.
Cao
, and
J. H.
Lu
,
SIAM J. Appl. Dyn. Syst.
7
,
108
(
2008
).
8.
J.
Zhou
and
J. A.
Lu
,
Physica A
386
,
481
(
2007
).
9.
K.
Park
,
Y.-C.
Lai
,
S.
Gupte
, and
J.-W.
Kim
,
Chaos
16
,
015105
(
2006
).
10.
W. L.
Lu
and
T. P.
Chen
,
Physica D
213
,
214
(
2006
).
11.
K. H.
Wang
,
X. C.
Fu
, and
K. Z.
Li
,
Chaos
19
,
023106
(
2009
).
12.
K. Z.
Li
,
M.
Small
,
K. H.
Wang
, and
X. C.
Fu
,
Phys. Rev. E
79
,
067201
(
2009
).
13.
K. Z.
Li
,
M. C.
Zhao
, and
X. C.
Fu
,
IEEE Trans. Circuits Syst., I: Regul. Pap.
56
,
2280
(
2009
).
14.
R.
Pastor-Satorras
and
A.
Vespignani
,
Phys. Rev. E
63
,
066117
(
2001
).
15.
M. E. J.
Newman
,
SIAM Rev.
45
,
167
(
2003
).
16.
M.
Small
,
D. M.
Walker
, and
C. K.
Tse
,
Phys. Rev. Lett.
99
,
188702
(
2007
).
17.
X.
Fu
,
M.
Small
,
D. M.
Walker
, and
H.
Zhang
,
Phys. Rev. E
77
,
036113
(
2008
).
18.
H. F.
Zhang
,
K. Z.
Li
,
X. C.
Fu
, and
B. H.
Wang
,
Chin. Phys. Lett.
26
,
068901
(
2009
).
19.
K. Z.
Li
,
M.
Small
,
H. F.
Zhang
, and
X. C.
Fu
,
Nonlinear Anal.: Real World Appl.
11
,
1017
(
2010
).
20.
S.
Funk
,
E.
Gilad
,
C.
Watkins
, and
V. A. A.
Jansen
,
Proc. Natl. Acad. Sci. U.S.A.
106
,
6872
(
2009
).
21.
B. T.
Grenfell
,
O. N.
Bjørnstad
, and
J.
Kappey
,
Nature
414
,
716
(
2001
).
22.
N. C.
Grassly
,
C.
Fraser
, and
G. P.
Garnett
,
Nature
433
,
417
(
2005
).
23.
G.
Yan
,
Z. Q.
Fu
,
J.
Ren
, and
W. X.
Wang
,
Phys. Rev. E
75
,
016108
(
2007
).
24.
D. J. D.
Earn
,
P.
Rohani
, and
B. T.
Grenfell
,
Proc. R. Soc. London, Ser. B
265
,
7
(
1998
).
25.
X.
Li
and
G.
Chen
,
IEEE Trans. Circuits Syst., I: Fundam. Theory Appl.
50
,
1381
(
2003
).
26.
M.
Krstic
,
I.
Kanellakopoulos
, and
P.
Kokotovic
,
Nonlinear and Adaptive Control Design
(
Wiley
,
New York
,
1995
).
27.
Z.
Li
,
L. C.
Jiao
, and
J. J.
Lee
,
Physica A
387
,
1369
(
2008
).
28.
A.-L.
Barabási
and
R.
Albert
,
Science
286
,
509
(
1999
).