There are certain correlations between collective behavior and spreading dynamics on some real complex networks. Based on the dynamical characteristics and traditional physical models, we construct several new bidirectional network models of spreading phenomena. By theoretical and numerical analysis of these models, we find that the collective behavior can inhibit spreading behavior, but, conversely, this spreading behavior can accelerate collective behavior. The spread threshold of spreading network is obtained by using the Lyapunov function method. The results show that an effective spreading control method is to enhance the individual awareness to collective behavior. Many real-world complex networks can be thought of in terms of both collective behavior and spreading dynamics and therefore to better understand and control such complex networks systems, our work may provide a basic framework.

From the real world, we can find many examples in which collective behavior and spreading behavior appear simultaneously and interplay with one another. In the stock market, the prices of different stocks will increase or decrease synchronously if relevant political canard spreads. Conversely, the price fluctuation also influences the spreading of relevant political canard. With the fast spread of an infectious disease in society, the times of avoiding assemblage and washing hand, etc., of people will increase in a synchronous way to protect themselves. On the other hand, this synchronous response will weaken the disease spread to some degree. In order to precisely control these collective and spreading behaviors and understand their interplay from the viewpoint of mathematics, the first step should be the construction of a suitable model which can display similar properties for these real dynamical networks. So, in this work, we will provide some mathematical models and address correlation between the collective and spreading dynamics on complex networks. The research results show that our models correspond closely with many real dynamical complex networks, and an effective spreading control method is to enhance the individual awareness to collective behavior.

It is well known that many real biological and social systems can be considered as dynamical complex networks. For examples, the daily activities of cows (eating, lying down, and standing) can be modeled by a dynamical network, where local dynamics of each cow is described by an oscillator of a piecewise linear dynamical system.1 The authors in Ref. 1 not only studied interesting dynamics such as synchronization but also developed some biological predictions. Many other examples can be found in Ref. 2. Now, let us ask: when an infectious disease spreads among these cows, what is the impact on their synchronous behavior? Can then this synchronous behavior weaken or strengthen the disease spread process? We feel that these problems can be resolved by examining the interplay between the collective and spreading dynamics on complex networks.

It is an interesting and important topic to consider the interplay between different dynamical behaviors appearing in complex networks. The correlation between traffic flow and epidemic spreading on complex networks was investigated numerically and theoretically in Ref. 3, for the case where the epidemic incidence was shaped by traffic-flow conditions and epidemic pathways were defined and driven by flows. The results in Ref. 3 provided a general framework for us to understand the spreading processes on complex traffic networks. We have investigated mathematically the correlation between the dynamical synchronization and the epidemic behavior on complex networks,4 and a very explicit condition for synchronization with respect to the epidemic rate was obtained. However, in this case, we only considered a special collective behavior, i.e., global synchronization. Moreover, the correlation between the dynamical synchronization and the epidemic behavior is unidirectional, i.e., the spreading behavior can influence the dynamical synchronization, but not vice versa. So, in this paper, we will extend our former research work to consider that the correlation between the collective behavior and the spread behavior is bidirectional and address further phase synchronization which may be a more general collective behavior.

In this work, we first propose a general bidirectional model of collective behavior and spreading dynamics on complex networks in Sec. II. Then in Secs. III–V, two concrete models are provided and studied, respectively, corresponding to two collective behaviors of dynamical behavior network, i.e., global synchronization and phase synchronization. By using the Lyapunov function method, we investigate the spread threshold of spreading on a network. In Sec. VI, we investigate finally the control problem of spreading behavior and provide an effective control strategy.

A general coupled model of collective and spreading behaviors on complex networks can be described as

Ẋ(t)=F(X(t),c(t)),Ẏ(t)=G(Y(t),E(t)),ċ(t)=H(Y(t),E(t)),
(1)

where X(t)=(x1(t),x2(t),,xN(t)) with xi(t)n denotes the state variable of the i-th individual in a dynamical behavior network with size N, which can exhibit collective behavior under suitable conditions, which is a necessary condition in this work. The coupling strength c(t)>0. The mapping F:(Nn,)Nn controls the dynamical change process of state variable X(t). In the second equality of Eq. (1), Y(t)d denotes the density variable of a spreading process on a network with maximal degree d. The variable E(t) is the error of collective behavior among state variables xi(t),i=1,2,,N and may be defined in different forms. The mapping G:(d,)d characterizes the dynamical change process of density variable Y(t). For the last equality, the function H:(d,) defines an adaptive law of the coupling strength c(t).

System (1) gives a bidirectional model between collective behavior and spread dynamics, where the dynamical behavior process X(t) can play a role in spread behavior Y(t) by embedding the error E(t), and the spread behavior Y(t) influences the dynamical behavior process X(t) by changing its coupling strength c(t). In real life, many dynamical phenomena can be described and explained by system (1). For example, when a political canard spreads, many relevant stocks will increase or decrease their prices synchronously for the explicit benefit of their corresponding corporations by closer communication. Conversely, the collective price fluctuation also accelerates or decelerates the spreading of relevant political canards. Similarly, when a certain infectious disease breaks out, people or animals will take some collective protective measures such as washing hands frequently, avoiding assemblage, resting frequently, etc. At the same time, these collective behaviors will further influence the disease spread. Fig. 1 gives a schematic diagram of three groups of collective and spreading behaviors which are interrelated closely, where the arrowhead shows the reliant relation between them.

FIG. 1.

Relationship between some groups of collective and spreading behaviors. For each group, bidirectional arrowhead implies interactional relationship. Many other probable cases are not listed and this is denoted by the ellipsis.

FIG. 1.

Relationship between some groups of collective and spreading behaviors. For each group, bidirectional arrowhead implies interactional relationship. Many other probable cases are not listed and this is denoted by the ellipsis.

Close modal

In Secs. III–V, we will consider two important collective behaviors, i.e., global synchronization and phase synchronization, and investigate the interplay between them and corresponding spreading behaviors. Undoubtedly, to deal with these problems the first step is that system (1) should be described precisely with the corresponding mappings F, G, H, and E.

Before constructing a concrete system, we should make the following basic assumptions. There is a weakly linear coupling between individuals in the dynamical behavior network in the beginning stage when spreading begins. And, there exists an interactional relationship between a dynamical behavior network and a spreading network, i.e., inhibiting or promoting each other.

Based on these assumptions, the model of SIS spread synchronization proposed in Ref. 4 and the framework of general model (1), we can construct a concrete system as

{ẋi(t)=f(xi(t))+c(t)j=1NaijHxj(t),İk(t)=λϕ(t)k[1Ik(t)]Θ(t)Ik(t),ċ(t)=βI(t)E1(t),
(2)

where i = 1, 2,…, N, k = 1, 2,…,d. Compared to the general model (1), we have X(t)=(x1(t),x2(t),,xN(t)) and Y(t)=(I1(t),I2(t),,Id(t))T, correspondingly. Moreover, xi(t)n denotes the state variable of the i-th node at the time t, and the function f(·) defines the local dynamics of each node and is supposed to be chaotic. The function c(t)>0 is the coupling strength and the matrix Hn×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,
(3)

where ki denotes the degree of node i. With these assumptions, the eigenvalues5 of matrix A can be given by 0=λ1λ2λN. The global synchronization error is set as

E(t)=1Nj=1Ns(t)xj(t)21+s(t)xj(t)2[0,1),
(4)

where s(t) is the synchronous state of the dynamical behavior network. Then, we define ϕ(t)=(1α)E(t)+α with constant α(0,1), and E1(t)=j=1Ns(t)xj(t)2 in the third equation of system (2).

The variables Ik(t) denote the density of infected nodes (individuals) with connectivity (contact) k at time t and I(t)=k=1dp(k)Ik(t) is the total infectious density. The spread rate λ(0,1] denotes the probability with which each susceptible node is infected if it is connected to one infected node. The term Θ(t) gives the probability that a randomly chosen link emanating from a node leads to an infected node. Moreover, Θ(t) has the form

Θ(t)=k=1dkp(k)Ik(t)k,
(5)

where the average degree k=k=1dkp(k). By this form, we mean that the connectivities of nodes in the spreading network are uncorrelated. The parameter β>0.

The initial condition of system (2) can be set as follows. The initial state xi(0) is chosen randomly from the real numbers and Ik(0)=ρk,c(0)=σ with 0<ρk1 and 0<σ1.

The physical meaning of model (2) was explained in detail in our former paper.4 Besides, the additional term ϕ(t)=(1α)E(t)+α in the second equation of model (2) denotes the admission rate,6 as the information of synchronization can be considered as a kind of individual awareness (or the risk perception). When all individuals achieve synchronization, i.e., E(t)0 as t, then the admission rate achieves the minimum α. Smaller value of parameter α means greater awareness to collective behavior. The case α=1 shows there is no awareness to the information of synchronization.

The infection control behavior of individuals within the community can be quantified by the variable ϕ(t)—this is the degree to which the behavior of the individual acts to reduce their risk of infection from others (and also risk of infecting others for a disease with a latent period): the rate of infection λ becomes ϕ(t)λ. Nonetheless, this individual behavior is a manifestation of the individuals behavior (through, for example, wearing of face masks, modification of hygiene practice, sharing of utensils, and washing of hands), and this is something which can be observed by others. There is a collectivization in ones response—if more people are wearing face masks in public it becomes more acceptable to do so and one is more likely to follow suit: or vice versa. Hence, the dynamical behavior of this parameter ϕ(t) exhibits a synchronization which, in term, influences the disease dynamics.

Now, we will address two basic properties of system (2), i.e., spread threshold of the spreading network and synchronization stability of the dynamical behavior network. By using a similar analytical method presented in Ref. 7, we can obtain easily that the spread threshold is

λc=kαk2>kk2.
(6)

In order to prove this spread threshold from the view of mathematics we should adopt a global stability analysis method, which will be shown in the next section. For the dynamical behavior network, if λ>λc, then all xi(t),i=1,2,,N will realize synchronization globally and asymptotically, as we will investigate in Sec. IV.

We now give some numerical examples to investigate model (2). Without loss of generality, we assume that the local dynamics of the dynamical behavior network are identical and defined as the chaotic Lorenz oscillation.8 From the point of view of a physical spread transmission process, it is rather difficult to justify this assumption. But we just use it for numerical simulations. This oscillation can be presented 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),
(7)

where parameters a1=10,a2=28,a3=(8/3). And the inner-coupling matrix H is chosen as the identity matrix.

First, we consider that the spreading network and dynamical behavior network having the same network topological structure. This implies that the individuals in the spreading network and dynamical behavior network are identical. For example, when human disease spreads, the reactive collective behaviors of washing hand, avoiding assemblage, etc., are also generated by people within the community. The network structure embedded in model (2) is set to be the BA network9 with size N = 200. This network evolves from initial network with size m0=4 and we add each new node with m = 3 new edges. Other parameters are chosen as ρ3=0.03,ρk=0,k3 and σ=0.01,β=0.001. Fig. 2 gives a simulation under λ=0.05,α=0.5, where the spreading does not become endemic and synchronization is not realized. In this case, the individuals do not exhibit collective behavior since the spread will die out and its prevalence is too small to keep them close. By increasing the spreading rate to λ=0.3 in Fig. 3, the synchronization is achieved and spreading also becomes endemic. Moreover, a very interesting phenomenon is that the total density I(t) is oscillating, which is also found in adaptive epidemic networks.10,11 If the collective behavior fails to influence the spreading, i.e., setting α=1, then the total density I(t) will converge monotonically to a larger value. From this simulation, we can see that the collective behavior can inhibit greatly the spreading behavior.

FIG. 2.

The changes of synchronization error E*(t)=1N1i=2Nx1(t)xi(t), epidemic prevalence I(t) and coupling strength c(t) in model (2) under parameters N=200,λ=0.05,α=0.5,β=0.001,σ=0.01. Both the spreading network and dynamical behavior network have the same BA network structure with minimal degree 3. There are only infected nodes with degree 3 and infection density ρ3=0.03, and the other nodes are all susceptible in the beginning stage.

FIG. 2.

The changes of synchronization error E*(t)=1N1i=2Nx1(t)xi(t), epidemic prevalence I(t) and coupling strength c(t) in model (2) under parameters N=200,λ=0.05,α=0.5,β=0.001,σ=0.01. Both the spreading network and dynamical behavior network have the same BA network structure with minimal degree 3. There are only infected nodes with degree 3 and infection density ρ3=0.03, and the other nodes are all susceptible in the beginning stage.

Close modal
FIG. 3.

The changes of E*(t), I(t) and c(t) in model (2) under parameters N=200,λ=0.3,α=0.5,β=0.001,σ=0.01.

FIG. 3.

The changes of E*(t), I(t) and c(t) in model (2) under parameters N=200,λ=0.3,α=0.5,β=0.001,σ=0.01.

Close modal

Next, we study an expert case in which the spreading network and dynamical behavior network have different network topological structures. This means that the individuals in the spreading network and the dynamical behavior network are nonidentical. For example, in the event of a political canard, the spreading network consists of individual people, while the stocks can be regarded as the corresponding dynamical behavior network of the stock market. Let N1,N2 denote the sizes of spreading network and dynamical behavior network both with BA mechanisms, respectively. With other parameters fixed, from Fig. 4 we can not find essential difference by only changing the network structure. The impact of the change of network structure on the spreading dynamics will be further discussed in detail in Sec. VI.

FIG. 4.

The changes of synchronization error E*(t)=1N21i=2N2x1(t)xi(t), epidemic prevalence I(t) and coupling strength c(t) in model (2) under parameters λ=0.3,α=0.5,β=0.001,σ=0.01. Both the spreading network and dynamical behavior network have different BA network structures with sizes N1=200,N2=300, respectively.

FIG. 4.

The changes of synchronization error E*(t)=1N21i=2N2x1(t)xi(t), epidemic prevalence I(t) and coupling strength c(t) in model (2) under parameters λ=0.3,α=0.5,β=0.001,σ=0.01. Both the spreading network and dynamical behavior network have different BA network structures with sizes N1=200,N2=300, respectively.

Close modal

In this section, we will study the global stability of equilibriums of model (2) by utilizing the method of global Lyapunov functions. Based on this analysis, we can get the spread threshold (6) for the spreading network and a synchronization condition for the dynamical behavior network.

For the dynamical behavior network in model (2), we first make the following preparations.

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

(xy)TP[f(x)f(y)]ξ(xy)T(xy).

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 error system of dynamical behavior network in model (2) can be written as

ė(t)=F¯(t)+c(t)(AH)e(t)+(INIn)G(t),
(8)

where is Kronecker product and IN denotes N-order identity matrix.

For the spreading network in model (2), we set βkj=λkjp(j)k,i,j=1,2,,d, and nonnegative matrices

M(I)=ϕ(t)(βkj(1Ik))d×d
(9)

and

M0=(βkj)d×d.
(10)

By setting I=(I1,I2,,Id)d, then the spreading network in model (2) can be rewritten as in a more compact form

İ(t)=M(I)II.
(11)

By simple computation, we know that the nonnegative matrix M0 has eigenvalues as μ1=μ2==μd1=0 and μd=λΣi=1di2p(i)k=λk2k. Define

R0=ϕ(t)ρ(M0),
(12)

where ρ denotes the spectral radius. Then we have R0=λϕ(t)k2k, which is the basic reproduction number for this spreading network that will be shown later. When the synchronization of the dynamical behavior network achieves stability, R0 will converges to fixed value λαk2k, from which we can get the spread threshold (6).

Obviously, the matrix M(I0) is irreducible, where I0=(0,0,,0)d. Define

Γ={(I1,I2,,Id)d|0Ik1,k=1,2,,d},

and let Γo denote the interior of Γ. The spreading network in model (2) is said to be uniformly persistent12 in Γo, if there exists a constant γ(0,1) such that liminftIk(t)>γ for all k provided (I1(0),I2(0),,Id(0))Γo. Since the spreading network can be reduced to a particular case of multi-group epidemic model, we have the following stability analysis mainly enlightened by the work in Ref. 12, which have solved the uniqueness and global stability of a multi-group SIR epidemic model.

Theorem 1. If R01, then I0 is the unique equilibrium of the spreading network in model (2) and it is globally stable in Γ. If R0>1, then I0 is unstable and this network is uniformly persistent in Γo.

Proof. By noting that

ρ(M(I))=ϕ(t)ρ{(βkj(1Ik))d×d}<ϕ(t)ρ(M0),

if R0=ϕ(t)ρ(M0)1, then the equation M(I) I = I has only the zero solution I=I0.

Let ω=(ω1,ω2,,ωd) be a left eigenvector of M0 corresponding to ρ(M0), i.e.,

ωρ(M0)=ωM0.
(13)

Since M0 is irreducible, then ωi>0 for i = 1, 2,…, d. Define the following function:

V(t)=k=1dωkIk.
(14)

The derivative of V(t) with respect to t along the solution of the system (11) is given by

dV(t)dt=ω[M(I)II]ω(ϕ(t)M0II)=ϕ(t)ωM0IωI=ϕ(t)ρ(M0)ωIωI=[ϕ(t)ρ(M0)1]ωI=(R01)ωI.
(15)

If R0<1, then dV(t)dt=0 means I=I0. If R0=1, then dV(t)dt=0 implies ωM(I)I=ωI. Assuming that II0, then we can get that ωM(I)I<ωϕ(t)M0I=ωI. So, in this situation, the equation ωM(I)I=ωI is satisfied if and only if I=I0. Therefore, when R01 the only compact invariant subset of the set {I|dV(t)dt=0} is the singleton {I0}. By LaSalle's Invariance Principle, I0 is globally asymptotically stable.

If R0>1 and II0, by using Eq. (15) we have

ω(ϕ(t)M0II)>0.

With the limitation limII0ω[M(I)II]=ω(ϕ(t)M0II), we can conclude that dV(t)dt>0 in a neighborhood of I0 in Γo. So, in this case the equilibrium I0 is unstable. By a similar discussion in Ref. 12, this instability means that this network is uniformly persistent in Γo.

Theorem 2. If R0>1, then there exists a unique endemic equilibrium I* of the spreading network in model (2), and it is globally asymptotically stable in Γo. Moreover, the synchronization manifold of the dynamical behavior network in this model is also globally asymptotically stable.

Proof. Since Sk(t)+Ik(t)=1 for k = 1, 2,…, d, the spreading network can be rewritten as

Ṡk(t)=1Skjϕ(t)βkjSkIj,
(16a)
İk(t)=Ik(t)+jϕ(t)βkjSkIj.
(16b)
And the error system of the dynamical behavior network is described by

ė(t)=F¯(t)+c(t)(AH)e(t)+(INIn)G(t).
(17)

Let (I1*,I2*,,Id*)Γo be an endemic equilibrium of the spreading network (16), and set β¯ij=βijSi*Ij*, where Si*(t)=1Ii*(t), then define a matrix

B¯=(j1β¯1jβ¯21β¯d1β¯12j2β¯2jβ¯d2β¯1dβ¯2djdβ¯dj),
(18)

whose each column sum equals zero.

Construct a function as V(t)=V1(t)+V2(t), with

V1(t)=k=1dvk(SkSk*lnSk+IkIk*lnIk),
(19)

where vk>0 is the cofactor of the k-th diagonal entry of B¯ satisfying B¯(v1,v2,,vd)TB¯v=0,12 and

V2(t)=12eT(t)(INP)e(t)+12β¯(c0c(t))2,
(20)

where β¯=λ2λmin(PH)>0,λmin(PH) denotes the minimal eigenvalue of matrix PH and c0 is a undetermined constant.

For V1(t), its derivative with respect to t along the solution of the system (16) is given by

dV1(t)dt=k=1dvk(1Skjϕ(t)βkjSkIjSk*Sk(1Skjϕ(t)βkjSkIj)Ik+jϕ(t)βkjSkIjjϕ(t)βkjSkIjIk*Ik+Ik*)=k=1dvk(Sk*+αjβkjSk*Ij*Skjϕ(t)βkjSkIjSk*Sk(Sk*+αjβkjSk*Ij*)+Sk*+jϕ(t)βkjSk*Ij+jϕ(t)βkjSkIjIkjϕ(t)βkjSkIjIk*Ik+αjβkjSk*Ij*)=k=1dvk(Sk*(SkSk*+Sk*Sk2)+2αjβkjSk*Ij*jϕ(t)βkjSkIjαjβkj(Sk*)2SkIj*+jϕ(t)βkjSkIjjϕ(t)βkjSkIjIk*Ik+jϕ(t)βkjSk*IjIk).
(21)

Now, we will show that

k=1dvk(j=1dϕ(t)βkjSk*IjIk)E1(t)Nk,j=1dvkβkjSk*Ij.
(22)

For the left-hand side of above inequality, we have

k=1dvk(j=1dϕ(t)βkjSk*IjIk)=k=1d(j=1dϕ(t)βjkSj*vjvk)IkE(t)k,j=1dvjβjkSj*Ik+k=1d(αj=1dβjkSj*vjvk)IkE1(t)Nk,j=1dvkβkjSk*Ij+k=1d(αj=1dβjkSj*vjvk)Ik.

To prove the inequality (22), it suffices to show

αj=1dβjkSj*vjvk=0,k=1,2,,d.

To this end, by using B¯v=0 we consider

αj=1dβjkSj*Ik*vj=αjkdβ¯jkvj+αβ¯kkvk=αjkdβ¯kjvk+αβ¯kkvk=αj=1dβ¯kjvk=αj=1dβkjSj*Ik*vk=Ik*vk,

that implies αj=1dβjkSj*vjvk=0 for all k. Since

SkSk*+Sk*Sk20,

we get

Sk*(SkSk*+Sk*Sk2)0,
(23)

and the equal sign holds if and only if Sk=Sk*.

Using Eqs. (21)–(23) and noting that ϕ(t)α>0 and β¯ij=βijSi*Ij*, we further obtain

dV1(t)dtE1(t)Nk,j=1dvkβkjSk*Ij+αk=1dvk(2j=1dβ¯kjj=1dβ¯kjSk*Skj=1dβ¯kjIjSkIk*IkSk*Ij*)=k,j=1dvkβkjSk*IjNi=1Nei(t)Tei(t)+αk,j=1dvkβ¯kj(2Sk*SkIjSkIk*IkSk*Ij*).
(24)

Based on graph theory, the authors in Ref. 12 have proven that k,j=1dvkβ¯kj(2Sk*SkIjSkIk*IkSk*Ij*)0 for positive βij. And the above equal sign holds if and only if Sk=Sk*,Ik=Ik*.

Now, let us turn to the derivative of V2(t) with respect to t along the solution of the system (16). By utilizing the similar analysis process present in Ref. 4, we can get

dV2(t)dt[ξ+c0λ2λmin(PH)βI(t)]i=1Nei(t)Tei(t).

From Theorem 1, we know that if R0>1, then the spreading network in model (2) is uniformly persistent in Γo. Combining the continuity and this uniformly persistent property of function I(t), we can conclude that if R0>1, then there is a constant γ¯>0 such that I(t)=k=1dp(k)Ik(t)>γ¯ for all t>0. So, we can further obtain

dV2(t)dt[ξ+c0λ2λmin(PH)βγ¯]i=1Nei(t)Tei(t).
(25)

Integrating the above discussions, we have

dV(t)dt[ξ+k,j=1dvkβkjSk*IjN+c0λ2λmin(PH)βγ¯]×i=1Nei(t)Tei(t).
(26)

Thus, we can select an adequately large constant c0 such that dV(t)dt0. Moreover, form inequality (24) and (26), we know the largest invariant subset of

{(S1,,Sd,I1,,Id,e1,,eN,c)|dV(t)dt=0}
(27)

is the singleton (S1*,,Sd*,I1*,,Id*,0,,0,c0). By LaSalle's Invariance Principle, this equilibrium is globally asymptotically stable. So, the unique endemic equilibrium I* of the spreading network in model (2) is globally asymptotically stable in Γo, and the synchronization manifold of the dynamical behavior network is also globally asymptotically stable.

In addition, note here that we have obtained the basic reproduction number by a global analysis to SIS model on complex network, while the results of the literature12 are applicable to multi-group SIR model. This means that we have extended the analysis in Ref. 12 to a more general model.

Compared to global synchronization discussed in Sec. IV, phase synchronization may be a more general collective behavior and more commonly observed in the real world. So, in this section, we will address this collective behavior and its influence on spreading behavior. Based on the famous Kuramoto model13–15 and the framework of general model (1), we can construct a concrete system as

{θ̇i(t)=ωi+c(t)j=1Naijsin(θj(t)θi(t)),İk(t)=λ[(1α)(1E(t))+α]k[1Ik(t)]Θ(t)          Ik(t),ċ(t)=βI(t)(1E(t)),
(28)

where i = 1, 2,…, N, k = 1, 2,…, d. The phase of the i-th individual is denoted by θi, and ωi represents its intrinsic frequency. Compared to the general model (1), we get X(t)=(θ1(t),θ2(t),,θN(t)) and Y(t)=(I1(t),I2(t),,Id(t))T, correspondingly. The phase synchronization error is set as

E(t)=|1Nj=1Neiθj|[0,1],
(29)

and the meaning of other mathematical symbols in model (28) is the same as that stated in Sec. III. If E(t)1 as t, then the dynamical behavior network achieves global phase synchronization. If E(t)0, then the phases of all individuals are different from each other and no synchronization phenomenon exists in this dynamical behavior network. When E(t)ε and 0<ε<1, this means that cluster synchronization appears with a proportion ε. Then, the spread threshold of the spreading network is

λc=1(1α)(1ε)+α×kk2>kk2.
(30)

The initial condition of system (28) can be set as follows. Without loss of generality, we choose the WS small-world network16 with probability p = 0.1 for rewiring links as the topology structure for the spreading network and dynamical behavior network. The initial phase θi(0) and intrinsic frequency ωi are chosen uniformly from intervals (π,π) and (−1/2, −1/2), respectively. The parameters α=0.8,β=1 and initial coupling strength c(0) = 0.1.

With the increasing spreading rates λ, we can see that the number of synchronous clusters in dynamical behavior network becomes smaller and the stable total density of spreading network becomes larger (see Figs. 5 and 6). Moreover, an interesting observation in Fig. 6 is that the total density fluctuates and goes to zero eventually. Then, by setting α = 1 we can see from Fig. 7 that the spreading behavior becomes endemic. From these simulations, we can conclude that the collective behavior of the dynamical behavior network can inhibit efficiently the spreading behavior. On the contrary, strong spreading behavior accelerates the collective behavior. These characteristics accord with many real dynamical networks very well.

FIG. 5.

The changes of synchronization error E(t), phase θi, epidemic prevalence I(t) and coupling strength c(t) in model (28) under parameters N=100,λ=0.05,α=0.8. Both the spreading network and dynamical behavior network have the same WS small-world network structure. There are only infected nodes with degree 3 and infection density ρ3=0.03, and the other nodes are all susceptible in the beginning stage.

FIG. 5.

The changes of synchronization error E(t), phase θi, epidemic prevalence I(t) and coupling strength c(t) in model (28) under parameters N=100,λ=0.05,α=0.8. Both the spreading network and dynamical behavior network have the same WS small-world network structure. There are only infected nodes with degree 3 and infection density ρ3=0.03, and the other nodes are all susceptible in the beginning stage.

Close modal
FIG. 6.

The changes of synchronization error E(t), phase θi, epidemic prevalence I(t) and coupling strength c(t) in model (28) under parameters N=100,λ=0.2,α=0.8.

FIG. 6.

The changes of synchronization error E(t), phase θi, epidemic prevalence I(t) and coupling strength c(t) in model (28) under parameters N=100,λ=0.2,α=0.8.

Close modal
FIG. 7.

The changes of synchronization error E(t), phase θi, epidemic prevalence I(t) and coupling strength c(t) in model (28) under parameters N=100,λ=0.2,α=1.

FIG. 7.

The changes of synchronization error E(t), phase θi, epidemic prevalence I(t) and coupling strength c(t) in model (28) under parameters N=100,λ=0.2,α=1.

Close modal

This section will address the control problem of the spreading network (2) by adjusting its structure and awareness to collective behavior and then provide an effective control strategy to prevent or weaken the diffusion of the spreading behavior. The results show that the awareness is a critical factor for this control strategy, while the network structure seems relatively insignificant in this control process.

The change of network structure is performed by adjusting the rewiring probability p in WS small-world network. By increasing the probability p from 0 to 1, we can get a transition from a regular network to a random graph. The awareness can be adjusted by changing the parameter α in model (2), where smaller value of parameter α means greater awareness to collective behavior.

As we known, the eigenvalue ratio ΛN/Λ2 of the adjacent matrix can quantify the synchronizability of the dynamical behavior network17,18 (the smaller this ratio is, the stronger synchronizability of the network), and the ratio k/k2 denotes the spread threshold of the traditional SIS network model (i.e., the second equation in model (2) with ϕ(t)=1). By increasing the probability p in WS small-world network, we find the ratios ΛN/Λ2 and k/k2 both decrease (see Fig. 8). This implies that the synchronizability is enhanced in the uncoupled dynamical behavior network, and the spread threshold decreases in the uncoupled spreading network. However, as these two networks are coupled by the form (2), we find that the change of network structure seems having relatively insignificant impact on the epidemics in this process, compared to the change of awareness.

FIG. 8.

Both eigenvalue ratio ΛN/Λ2 and spread threshold k/k2 decrease with increasing rewiring probability p in WS small-world network. This figure is obtained by 50 independent realizations.

FIG. 8.

Both eigenvalue ratio ΛN/Λ2 and spread threshold k/k2 decrease with increasing rewiring probability p in WS small-world network. This figure is obtained by 50 independent realizations.

Close modal

Initial condition setting: The network size N = 100, α=0.5,σ=0.01,β=0.001. The initial infectious density Idmin(0)=1Np(dmin),Ik(0)=0 for all kdmin (i.e., just an individual with minimal degree dmin is infected). All simulations in this section are based on 50 independent realizations.

From Figs. 9 and 10, we can see that under fixed parameter α = 0.5, the change of network structure just plays trivial role in the spreading process, as the spreading prevalence I(t) does not vary obviously either in disease-free case or in endemic case (By performing simulations with more nodes, we have not found any essential difference to this simulation with 100 nodes. So, we can conclude the qualitative dependence of I(t) on p by this two figures). For other values of parameter α, the results are similar to this case. We can further observe this trivial influence in Fig. 11, where the spread threshold λc seems fixed if the parameter α keeps constant. However, with the increasing of parameter α, the spread threshold λc is decreased greatly. Therefore, the awareness is a critical factor for the spreading control strategy, and an effective control method is to enhance the awareness to collective behavior.

FIG. 9.

The epidemic prevalence of model (2) with disease-free equilibrium under λ = 0.2. The rewiring probability p increases form 0 to 0.45 with step length 0.05. The color point denotes the value of total infectious density I(t).

FIG. 9.

The epidemic prevalence of model (2) with disease-free equilibrium under λ = 0.2. The rewiring probability p increases form 0 to 0.45 with step length 0.05. The color point denotes the value of total infectious density I(t).

Close modal
FIG. 10.

The epidemic prevalence of model (2) with endemic equilibrium under λ = 0.4. The rewiring probability p increases form 0 to 0.45 with step length 0.05.

FIG. 10.

The epidemic prevalence of model (2) with endemic equilibrium under λ = 0.4. The rewiring probability p increases form 0 to 0.45 with step length 0.05.

Close modal
FIG. 11.

The spread threshold under different awareness degree α and rewiring probability p.

FIG. 11.

The spread threshold under different awareness degree α and rewiring probability p.

Close modal

In summary, this paper has constructed several coupled models which can simulate collective and spreading behaviors on complex networks. Two concrete models are studied, respectively, where the spreading behavior is controlled by traditional SIS network model and the collective behavior is demonstrated by global synchronization and phase synchronization. The spread threshold of spreading network is obtained by using the stability theory method, and it depends on the network structure and individual awareness. The synchronization manifold of the dynamical behavior network is globally asymptotically stable if the spreading network can achieve an endemic state. Moreover, some numerical simulations are given to verify these theoretical results. We find that the collective behavior can inhibit spreading behavior, but, conversely, this spreading behavior can accelerate collective behavior. Finally, we study the impact of the change of network structure on spreading dynamics and find an effective control of spreading behavior is to enhance the awareness to collective behavior. This work may provide a basic framework to better understand and control such complex networks systems.

This research was supported jointly by NSFC grants (Nos. 61004101, 60964006, 11061010, 11072136, 11161013, and 61164020). The authors gratefully acknowledge the support of the Guangxi Natural Science Foundation Program (Nos. 2011GXNSFB018059 and 2011GXNSFA018136), and Shanghai University Leading Academic Discipline Project (A.13-0101-12-004).

1.
J.
Sun
,
E. M.
Bollt
,
M. A.
Porter
, and
M. S.
Dawkins
, “
A mathematical model for the dynamics and synchronization of cows
,”
Physica D
240
,
1497
(
2011
).
2.
M. E. J.
Newman
, “
The structure and function of complex networks
,”
SIAM Rev.
45
,
167
(
2003
).
3.
S.
Meloni
,
A.
Arenas
, and
Y.
Moreno
, “
Traffic-driven epidemic spreading in finite-size scale-free networks
,”
Proc. Natl. Acad. Sci. U.S.A.
106
,
16897
(
2009
).
4.
K. Z.
Li
,
X. C.
Fu
,
M.
Small
, and
Z. J.
Ma
, “
Adaptive mechanism between dynamical synchronization and epidemic behavior on complex networks
,”
Chaos
21
,
033111
(
2011
).
5.
X. F.
Wang
and
G. R.
Chen
, “
Pinning control of scale-free dynamical networks
,”
Physica A
310
,
521
(
2002
).
6.
Q. C.
Wu
,
X. C.
Fu
,
M.
Small
, and
X. J.
Xu
, “
The impact of awareness on epidemic spreading in networks
,”
Chaos
22
,
013101
(
2012
).
7.
X. C.
Fu
,
M.
Small
,
D. M.
Walker
, and
H. F.
Zhang
, “
Epidemic dynamics on scale-free networks with piecewise linear infectivity and immunization
,”
Phys. Rev. E
77
,
036113
(
2008
).
8.
D. M.
Li
,
J. A.
Lu
,
X. Q.
Wu
, and
G. R.
Chen
, “
Estimating the ultimate bound and positively invariant set for the Lorenz system and a unified chaotic system
,”
J. Math. Anal. Appl.
323
,
844
(
2006
).
9.
A.-L.
Barabási
and
R.
Albert
, “
Emergence of scaling in random networks
,”
Science
286
,
509
(
1999
).
10.
T.
Gross
,
C.
Dommar D'Lima
, and
B.
Blasius
, “
Epidemic dynamics on an adaptive network
,”
Phys. Rev. Lett.
96
,
208701
(
2006
).
11.
L. B.
Shaw
and
I. B.
Schwartz
, “
Fluctuating epidemics on adaptive networks
,”
Phys. Rev. E
77
,
066101
(
2008
).
12.
H. B.
Guo
,
M. Y.
Li
, and
Z. S.
Shuai
, “
Global stability of the endemic equilibrium of multigroup SIR epidemic models
,”
Can. Appl. Math. Q.
14
,
259
(
2006
).
13.
S. H.
Strogatz
,
R. E.
Mirollo
, and
P. C.
Matthews
, “
Coupled nonlinear oscillators below the synchronization threshold: Relaxation be generalized landau damping
,”
Phys. Rev. Lett.
68
,
2730
(
1992
).
14.
H.
Hong
,
M. Y.
Choi
, and
B. J.
Kim
, “
Synchronization on small-world networks
,”
Phys. Rev. E
65
,
026139
(
2002
).
15.
M. G.
Earl
and
S. H.
Strogatz
, “
Synchronization in oscillator networks with delayed coupling: A stability criterion
,”
Phys. Rev. E
67
,
036204
(
2003
).
16.
D. J.
Watts
and
S. H.
Strogatz
, “
Collective dynamics of small world networks
,”
Nature
393
,
440
(
1998
).
17.
L.
Donetti
,
P. I.
Hurtado
, and
M. A.
Muñoz
, “
Entangled networks, synchronization, and optimal network topology
,”
Phys. Rev. Lett.
95
,
188701
(
2005
).
18.
C. Y.
Yin
,
W. X.
Wang
,
G. R.
Chen
, and
B. H.
Wang
, “
Decoupling process for better synchronizability on scale-free networks
,”
Phys. Rev. E
74
,
047102
(
2006
).