- Email: [email protected]
Epidemic spreading over quenched networks with local behavioral response
Chaos, Solitons and Fractals 96 (2017) 17–22
Contents lists available at ScienceDirect
Chaos, Solitons and Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos
Epidemic spreading over quenched networks with local behavioral response Qingchu Wu a,∗, Shufang Chen b, Lingling Zha a a b
College of Mathematics and Information Science, Jiangxi Normal University, Nanchang, Jiangxi 330022, P. R. China College of Physics and Communication Electronics, Jiangxi Normal University, Nanchang, Jiangxi 330022, P. R. China
a r t i c l e
i n f o
Article history: Received 14 October 2016 Revised 20 December 2016 Accepted 8 January 2017
Keywords: Complex networks Epidemic spreading Pair approximation
a b s t r a c t We discuss the impact of local information-based behavioral response on epidemic spreading in social networks. By using a pair quenched mean-field approach developed by Mata and Ferreira [Europhys. Lett. 103 (2013) 48003], we derive a dynamical model governing the epidemic spreading over a random network with a linear response function and density-dependent epidemic information. A deterministic relation between the epidemic threshold and the response parameter is derived according to a quasi-static approximation method. It is found that local behavioral response will induce the extinction of the disease via rasing the epidemic threshold. Additionally, the theoretical result is supported by stochastic simulations on an Erdo¨ s–Rényi random network and a Baraba´ si–Albert scale-free network. Simulations show that the pair quenched mean-field approach is more accurate than the classical quenched mean-field approach. © 2017 Elsevier Ltd. All rights reserved.
1. Introduction The theory of complex network involves the network itself (structure and its dynamics), dynamics on networks and coupled interactions between network and dynamics evolving over the network [1]. As a typical process occurring on the network [2], disease transmission dynamics have drawn a wide attention of researchers from mathematical, computer, physical and biological communities [3–6]. Recently, many researches focus on the coupled disease-behavior dynamics [7]. Generally speaking, individual behavioral response means that individual adjusts its adaptive behavior by changing contact numbers or contact objects [8,9], vaccination decision-Erdo¨ s-Rényimaking [10,11] and dynamical parameters [12–14]. In reality, the epidemic information (denoted by x) can cause individuals to keep social distancing (by wearing protective masks, washing hands frequently, avoiding crowded public areas, or staying home from work or school or more creative precautions), which potentially results in the reduction of individual susceptibility (denoted by y) [13,16]. For a node with k contact number and s infected neighbors, if the epidemic information obtained by a node is determined by the infection fraction in its neighborhood, then
∗
Corresponding author: E-mail address: [email protected] (Q. Wu).
http://dx.doi.org/10.1016/j.chaos.2017.01.003 0960-0779/© 2017 Elsevier Ltd. All rights reserved.
x = s/k and we say it is density-dependent [15–18]; if the information is only related to the number of infected neighbors, then x = s and it is frequency-dependent [19]. To incorporate the information effect into the transmission model, it is needed to define a response function y = φ (x ), which satisfies φ (0 ) = 0, 0 ≤ φ (x ) ≤ 1, and φ (x) < 0. It is relevant to investigate the effect of individual behavioral response on the epidemic spreading under the network framework with the spatial diffusion [20–22]. In [16], we studied the susceptible-infected-susceptible (SIS) model with a response function of linear form, i.e., φ (x ) = 1 − α x(1 ≥ α ≥ 0 ). Meanwhile, the epidemic information is densitydependent. By using the heterogeneous mean-field (HMF) approach, we obtained the epidemic threshold above which the epidemic spreads over the whole network, depending on the response parameter α . Shang [17] analyzed the case of nonlinear form φ (x ) = 1 − α (s/k )m (m is a positive integer) and derived the epidemic threshold. Recently, Zhang et al. [19] proposed a response function of exponential form φ (x ) = (1 − α )x with the frequency-dependent epidemic information. Their theoretical analysis and simulations show that the epidemic threshold is strongly affected by individual response. So, the spread of an infectious disease over a network can be controlled by individual behaviorial response. Compared to the immunization and quarantine, individual response is a very economic strategy. All these researches are based on the HMF model, that has been shown to be exact only for the annealed network [23,24]. How to analyze this issue on the static or quenched network? [25] How to
18
Q. Wu et al. / Chaos, Solitons and Fractals 96 (2017) 17–22
study the interplay between the network structure and the behavioral response? For the SIS-type model, the microscopic Markov-chain (MMC) approach is frequently-used in the analysis of the epidemic spreading on the quenched network [26]. Van Mieghem et al. developed a continuous-time MMC approach and called it as the N-intertwined model [5], which is also referred to as the quenched mean-field (QMF) model [28]. Compared to the HMF modeling framework, the N-interwind model may be more adapted for the quenched network [24]. However, this model is not enough exact since it only applies one-order neighbor information. Therefore, a highorder mathematical formulation is required to provide a more exact threshold analysis [27]. Mata and Ferreira proposed a pair-type formulation of quenched mean-field model, referred to as the pair quenched mean-field model [28]. The pair QMF model makes use of the two-order neighbor information and gives a more exact prediction. More importantly, the epidemic threshold can be determined and analyzed by this model. In the present work, we would like to give a rigorous derivation for the pair QMF model, and show that the pair QMF model is adapted for the static randomly connected network with no clustering [28]. In addition, we analyze the impact of local behavioral response on the epidemic spreading over the static network by using the pair QMF model, which is still challenging for us since it involves the relation between the contact structure and local behavioral response.
Within this perspective, they built a system describing the SIS epidemic dynamics on the network G, which takes the form
2. An analysis framework
2.2. The derivation of the pair QMF equations
We assume that an infectious disease follows SIS dynamics on a given network with a size N, denoted by G = (V, E ). Generally, the contact network is weighted [29] or directed [30]. For simplicity, we assume that G is a unweighted and undirected network and hence it is completely determined by its adjacency matrix A = (ai j ): if node i links to node j in G, then ai j = 1; otherwise ai j = 0. In this model, each node may stay in either susceptible (S) or infected (I) state. During a infinitesimal time interval (t , t + t ], an infected node transmits the infectious disease to its susceptible neighbor with probability β t and meanwhile, it recovers and becomes susceptible again with probability γ t. As usual, we define λ = β /γ be the effective spreading rate [4,25]. The mathematical model for the static contact network can be built by the microscopic Markov-chain approach [26,31]. This approach focuses on the dynamic of the probability of each node i to be infected at time t, denoted by ρ i . Let Xi (t) ∈ {0, 1} denote the state of node i at time t [5]. If Xi (t ) = 0, node i is susceptible; if Xi (t ) = 1, node i is infectious. Then ρi = P[Xi (t ) = 1] and 1 − ρi = P[Xi (t ) = 0], where P[A] represents the probability of event A occurring.
We will apply the total probability formula to derive the pair QMF equations for the SIS epidemic model (1). For the sake of the following analysis, we first present a lemma.
2.1. The pair quenched mean-field model Mata and Ferreira [28] proposed a pair-type formulation of quenched mean-field model, i.e., the pair QMF model. In order to develop the pair QMF theory, they introduced a series of simple notations: [Ai ] = P[Xi (t ) = A], which is the probability that node i is in the state A; [Ai , B j ] = P[Xi (t ) = A, X j (t ) = B], which is the probability that nodes i and j are in states A and B, respectively; [Ai , Bj , Cl ] is the extension to three nodes. In addition, they also defined specific pair-type variables: ωi j = [0i , 0 j ], φi j = [0i , 1 j ], φ¯ i j = [1i , 0 j ], ψi j = [1i , 1 j ]. It is easy to know that pair-type variables must satisfy
ψi j = ρ j − φi j , φ¯ i j = ρi − ψi j = ρi − ρ j + φi j , ωi j = 1 − ρi − φi j .
d ρi = −γ ρi + β φi j a i j dt N
j=1
d φi j = −γ φi j − βφi j + γ ψi j dt +β [0i , 0 j , 1l ] − β [1l , 0i , 1 j ] l∈N ( j ) l=i
(1)
l∈N (i ) l= j
where N (i ) denotes the neighborhood of node i. In order to close the above equations, the authors used the standard approximation [28,32]
[Ai , B j , Cl ]
[Ai , B j ][B j , Cl ] . [B j ]
(2)
In Appendix A, we prove that this approximation holds for the random network with no clustering. After performing a quasi-static approximation for t → ∞, the authors found that the epidemic threshold can be determined by the largest eigenvalue of the Jacobian matrix L, which is given by
Li j = −
γ+
β 2 ki β ( 2γ + β ) ai j . δi j + 2β + 2γ 2β + 2γ
(3)
Here, δ ij is the Kronecker symbol and ki is the degree of node i.
Lemma 1. For each a subset of nodes in G, V∗ (i.e., V∗ ⊂ V), we assume that a probability of each node i ∈ V∗ to be infected, denoted by ρ i is given, then the number of infected nodes in this subset is a stochastic variable ξ ∈ [0, |V∗ |] (|V∗ | denotes the number of the elements in set V∗ ) and its first and second order moment satisfy
E[ξ ] =
ρi , E [ ξ 2 ] =
i∈V ∗
ρi +
i∈V ∗
ρi 1 ρi 2 ,
i1 =i2
respectively. This lemma is easily proved by the mathematical induction. In fact, the first order moment formula has been proved in [31,33]. The second order formula can be justified by a simple example when |V ∗ | = 3. At this time, E[ξ 2 ] = 12 × [ρ1 (1 − ρ2 )(1 − ρ3 ) + ρ2 (1−ρ1 )(1 − 3 ) + ρ3 (1−ρ1 )(1−ρ2 )] + 22 × [ρ1 ρ2 (1 − ρ3 ) + ρ1 ρ3 (1 − ρ2 ) + ρ2 ρ3 (1 − ρ1 )] + 32 × ρ1 ρ2 ρ3 = ρ1 + ρ2 + ρ3 + 2(ρ1 ρ2 + ρ2 ρ3 + ρ1 ρ3 ), which accords with the conclusion from Lemma 1. In the following, we derive system (1). By the total probability formula, we have
P[Xi (t + t ) = 1] = P[Xi (t ) = 1]P[Xi (t + t ) = 1|Xi (t ) = 1] + P[Xi (t ) = 0]P[Xi (t + t ) = 1|Xi (t ) = 0].
(4)
Denote kinf (i) be the number of infected nodes in the neighborhood of node i. Note that
P[Xi (t ) = 0]P[Xi (t + t ) = 1|Xi (t ) = 0] = P[Xi (t + t ) = 1, Xi (t ) = 0] =
ki
P[Xi (t + t ) = 1, Xi (t ) = 0, kinf (i ) = s]
s=0
=
ki s=0
P[Xi (t + t ) = 1|Xi (t ) = 0, kinf (i ) = s]
Q. Wu et al. / Chaos, Solitons and Fractals 96 (2017) 17–22
P[Xi (t + t ) = 0, X j (t + t ) = 1|Xi (t ) = 1, X j (t ) = 0] = o(t ).
P[Xi (t ) = 0, kinf (i ) = s] ki
=
[1 − (1 − β t )s ]P[Xi (t ) = 0, kinf (i ) = s]
s=0
= E pa [1 − (1 − β t )s ] Then, we have
P[Xi (t + t ) = 1] = P[Xi (t ) = 1]P[Xi (t + t ) = 1|Xi (t ) = 1] +
ki
[1 − (1 − β t )s ]P[Xi (t ) = 0, kinf (i ) = s]
s=0
= P[Xi (t ) = 1](1 − γ ) +
ki
[1 − (1 − β t )s ]P[Xi (t ) = 0, kinf (i ) = s]
s=0
Thus, subtracting P[Xi (t ) = 1] and dividing by t and letting t → 0, we get
d ρi = −γ ρi + β dt
ki
sP[Xi (t ) = 0, kinf (i ) = s].
s=0
Xi (t ) = 0, X j (t ) = 0, kinf ( j ) = s]
β t P[Xi (t ) = 0, X j (t ) = 1, kinf (i ) − 1 = s]
ki −1
(5)
P[Xi (t + t ) = 1, X j (t + t ) = 1, Xi (t ) = 0, X j (t ) = 1,
s=0
kinf (i ) = s] ki −1
=
P[Xi (t + t ) = 1, X j (t + t ) = 1|Xi (t ) = 0, X j (t ) = 1,
s=0
kinf (i ) = s] P[Xi (t ) = 0, X j (t ) = 1, kinf (i ) = s] ki −1
=
[1 − (1 − β t )s ]P[Xi (t ) = 0, X j (t ) = 1, kinf (i ) = s]
s=0
k j −1
P[Xi (t + t ) = 0, X j (t + t ) = 1|Xi (t ) = 0,
s=0
X j (t ) = 0, kinf ( j ) = s] P[Xi (t ) = 0, X j (t ) = 0, kinf ( j ) = s] k j −1
[1 − (1 − β t )s ]P[Xi (t ) = 0, X j (t ) = 0, kinf ( j ) = s]
s=0
s = Etr 1 [1 − (1 − β t ) ],
and
P[Xi (t + t ) = 0, X j (t + t ) = 1|Xi (t ) = 1, X j (t ) = 1] = γ t + o(t ), and
ki −1
=
It is needed to give the specific expression of the four terms on the right-hand-side (RHS) of Eq. (5). For the former three terms, we have
P[Xi (t + t ) = 0, X j (t + t ) = 1,
P[Xi (t ) = 0, X j (t ) = 1, kinf (i ) = s]
P[Xi (t + t ) = 1, X j (t + t ) = 1, Xi (t ) = 0, X j (t ) = 1]
+ P[Xi (t + t ) = 0, X j (t + t ) = 1|Xi (t ) = 0, X j (t ) = 1]
k j −1
kinf (i ) − 1 = s]
For part (ii), we turn to compute the probability P[Xi (t + t ) = 1, X j (t + t ) = 1, Xi (t ) = 0, X j (t ) = 1], which obeys
P[Xi (t ) = 1, X j (t ) = 1]
P[Xi (t + t ) = 0, X j (t + t ) = 1, Xi (t ) = 0, X j (t ) = 0]
P[Xi (t + t ) = 1, X j (t + t ) = 1|Xi (t ) = 0, X j (t ) = 1,
s=0
+ P[Xi (t + t ) = 0, X j (t + t ) = 1|Xi (t ) = 1, X j (t ) = 1]
P[Xi (t ) = 0, X j (t ) = 1]
= Etr 2 [β t] = β t P[Xi (t ) = 0, X j (t ) = 1].
+ P[Xi (t + t ) = 0, X j (t + t ) = 1|Xi (t ) = 1, X j (t ) = 0]
P[Xi (t + t ) = 1, X j (t + t ) = 1, Xi (t ) = 0, X j (t ) = 1,
s=0
P[Xi (t ) = 1, X j (t ) = 0]
=
ki −1
=
=
= P[Xi (t + t ) = 0, X j (t + t ) = 1, Xi (t ) = 0, X j (t ) = 0]
ki −1
kinf (i ) − 1 = s]
s=0
P[Xi (t + t ) = 0, X j (t + t ) = 1]
=
P[Xi (t + t ) = 1, X j (t + t ) = 1, Xi (t ) = 0, X j (t ) = 1]
s=0
Following Lemma 1, s=0 sP[Xi (t ) = 0, kinf (i ) = s] = j∈N (i ) P[Xi (t ) = 0, X j (t ) = 1]. So the first equation of Eq. (1) is derived. From this equation, one finds that the term φ ij has no influence on the dynamic of ρ i when node i does not connect to node j. Hence, it is required for us to further analyze the dynamic of φ ij as a pair of connected node (i, j). By the total probability formula, we have
While, the computation of the fourth term in the RHS of Eq. (5) is not straightforward. We notice that P[Xi (t + t ) = 0, X j (t + t ) = 1|Xi (t ) = 0, X j (t ) = 1] = 1−P[Xi (t + t ) = 0, X j (t + t ) = 0|Xi (t ) = 0, X j (t ) = 1] − P[Xi (t + t ) = 1, X j (t + t ) = 0|Xi (t ) = 0, X j (t ) = 1] − P[Xi (t + t ) = 1, X j (t + t ) = 1|Xi (t ) = 0, X j (t )=1]. Clearly, P[Xi (t + t ) = 0, X j (t + t ) = 0|Xi (t ) = 0, X j (t ) = 1] = γ t and P[Xi (t + t ) = 1, X j (t + t ) = 0|Xi (t ) = 0, X j (t ) = 1] = o(t ). It is a little complicated for the computation of P[Xi (t + t ) = 1, X j (t + t ) = 1|Xi (t ) = 0, X j (t ) = 1] because we do not know about the number of infected neighbors of node i. To solve this, we decompose it to two parts: (i) node i is infected by node j since j is infectious; (ii) node i is infected by node l (l = j. For part (i), it is easy to get that P[Xi (t + t ) = 1, X j (t + t ) = 1|Xi (t ) = 0, X j (t ) = 1] = β t. However, this also can be obtained by using the property of the joint probability:
=
k i
=
19
s = Etr 3 [1 − (1 − β t ) ].
Finally, subtracting P[Xi (t ) = 0, X j (t ) = 1] and dividing by t and letting t → 0 in Eq. (5) and using Lemma 1, one can obtain the second equation of Eq. (1). It is worthy noting that the derivation of Eq. (1) does not need any external assumption, so it is true for any network structure. 3. The local behavioral response model 3.1. The model The objective of the present work is to investigate the effect of the local epidemic information on the epidemic dynamic of SIS model. In our model, it is assumed that the response function is a linear formulation with respect to the epidemic information and is also density-dependent [16]. It suffices to replace β
20
Q. Wu et al. / Chaos, Solitons and Fractals 96 (2017) 17–22
by β [1 − α kinf (i )/ki ] for susceptible node i. We apply the analysis framework as mentioned in Section 2 to derive corresponding pair QMF equations. Note that
E [1 − (1 − β (1 − α s/ki )t ) ] t pa
lim
t→0
s
βα
= β E pa [s] −
E pa [s2 ] ki α βα =β 1− [0i , 1 j ] − ki ki j∈N (i )
N α
=β 1− lim
t→0
ki
φi j a i j − β
j=1
ki
=β
=β 1−
−
α
βα
[0i , 1 j1 ][0i , 1 j2 ]
d α ρi = −γ ρi + β 1 − dt ki
j1 , j2 ∈N (i ), j1 = j2
φi j 1 φi j 2 a i j 1 a i j 2 ,
kj
+β 1 −
[0i , 0 j , 1l ] −
l∈N ( j ) l=i
βα kj
[0i , 0 j , 1l1 ][0i , 0 j , 1l2 ],
l1 ,l2 ∈N ( j ) l1 ,l2 =i,l1 =l2
Etr βα tr 2 [β (1 − α (s + 1 )/ki )t] = β Etr E [s] 2 [1] − t→0 t ki 2 α βα =β 1− [1 , 0i , 1 j ] φi j − ki ki l∈N (i) l lim
l= j
and
lim
t→0
s Etr 3 [1 − (1 − β (1 − α s/ki )t ) ] t
= β Etr 3 [s] −
Etr [s2 ] ki 3 α βα =β 1− [1 , 0i , 1 j ] − ki l∈N (i) l ki l= j
d φi j = dt
β 1− −
βα kj
α kj
N
[1l1 , 0i , 1 j ][1l2 , 0i , 1 j ].
l1 ,l2 ∈N (i ) l1 ,l2 = j,l1 =l2
φi j a i j − β
j=1
+
βα ki
ki
φi j 1 φi j 2 a i j 1 a i j 2
ki
j
j
2γ + β (1 − kα ) + β (1 − kα )
.
(8)
j
ui j =
β (1 − α /k j ) , 2γ + β (1 − α /ki ) + β (1 − α /k j )
vi j =
2γ + β (1 − α /k j ) . 2γ + β (1 − α /ki ) + β (1 − α /k j )
After plugging (8) in (7), we obtain
d ρi = L˜i j ρ j , dt
j1 = j2
where the system Jacobian is given by
⎛ βα ⎝ ki
1+
[1l , 0i , 1 j ]⎠
l∈N (i ) l= j
[1l , 0i , 1 j ]
l∈N (i ) l= j
[1l1 , 0i , 1 j ][1l2 , 0i , 1 j ]
l1 ,l2 ∈N (i ) l1 ,l2 = j,l1 =l2
ωi j φ jl φ¯ li φi j , [1l , 0i , 1 j ] = . 1 − ρj 1 − ρi
γ +β 1−
N α
ki
ai j ui j
α δi j + β 1 − ai j vi j . ki
j=1
˜ PQMF , is deterAs a result, the epidemic threshold, denoted by λ c mined by the equation max ( L ) = 0. Here, max (M) represents the largest eigenvalue of the matrix M. So, the outbreak condition can be seen to depend only of the values of α , A and β /γ . Finally, we analyzed the local response model by taking an one-vertex approximation, i.e., P[Xi (t ) = 0, X j (t ) = 1] = P[Xi (t ) = 0]P[X j (t ) = 1], or φi j (1 − ρi )ρ j . At this time, the epidemic spreading is simply governed by
⎞
(9)
j=1
[0i , 0 j , 1l1 ][0i , 0 j , 1l2 ]
α
[2γ + β (1 − kα )]ρ j − β (1 − kα )ρi
Li j = −
where
[0i , 0 j , 1l ] =
α
[0i , 0 j , 1l ]
+ γ ψi j − γ φi j − βφi j +
(7)
N
l1 ,l2 ∈N ( j ) l1 ,l2 =i,l1 =l2
−β 1 −
φ jl (a jl − δil )
l=1
For clarity, let
l∈N ( j ) l=i
kj
N
φi j + γ ρ j
j
Upon substituting these into the Eqs. (4) and (5), we write down equations for the changes in the SIS model with local behavioral response
d α ρi = −γ ρi + β 1 − dt ki
α
φi j a i j
j=1
However, it seems hard to derive the epidemic threshold for system (7) according to the stability of its Jacobin matrix [34], due to its complicated formulation. Indeed, as for the traditional pairapproximation model on heterogeneous networks [35,36], it is not easy to explicitly derive analytical expressions of the outbreak condition. However, following [28,32,37], a quasi-static approximation can be applied into our model to solve this problem. For t → ∞, we let dρ i /dt ࣃ 0 and dφ i, j /dt ࣃ 0 when the system is around the zero solution. Hence, we have β (1 − kα ) j N N l=1 φ jl a jl = γ ρ j and l=1 φ jl δil = φ ji . Then, substituting these equations in Eq. (7) leads to
φi j =
βα
N
d α φi j = − 2 γ + β 1 − dt ki
j1 = j2
2 Etr 1 [s ]
kj
j
ψi j + φi j = ρ j , we have
s Etr 1 [1 − (1 − β (1 − α s/k j )t ) ] t
Etr 1 [s]
In the epidemiology, it is critical to determine the epidemic threshold, by which one can judge whether the infectious disease spreads in a population. To this end, we write a linear form of the model around the fixed point ρi = φi j = ψi j = 0. At this time, N ωi j φ jl N l∈N ( j ) [0i , 0 j , 1l ] = l=1 1−ρ (a jl − δil ) l=1 φ jl (a jl − δil ). Since l=i
α
3.2. Epidemic threshold
(6)
d ρi = −γ ρi dt + β ( 1 − ρi )
1−
N α
ki
j=1
ρ j ai j −
α ki
ρ j1 ρ j2 ai j1 ai j2 .
j1 = j2
(10) Denote Li j be the Jacobian matrix of Eq. (10). When its largest ˆ QMF ) eigenvalue is null, the epidemic threshold (denoted by λ c
Q. Wu et al. / Chaos, Solitons and Fractals 96 (2017) 17–22
21
L + 1 and final infection size ρ as functions of the infection rate β . Rpqmf = 1 and Rqmf = 1 imply that max (L˜ ) = 0 Fig. 1. Illustrations of Rpqmf = max ( L ) + 1, Rqmf = max L = 0, respectively. They represent the critical levels of epidemic outbreak. and max
satisfies
λˆ QMF = c
1
max (H )
,
(11)
where H = [hi j ]N×N , hi j = ai j − α ai j /ki . 3.3. Simulations To test above argument, we perform stochastic simulations over both two typical contact networks. One is a Baraba´ si–Albert (BA) scale-free network with the average degree k = 6 [39], which has a highly heterogeneous degree distribution (power-law distribution) and a small clustering coefficient. The other is an Erdo¨ s–Rényi (ER) random network with connecting probability p = 0.006 [38], which has a homogeneous degree distribution (Poisson distribution) and a small clustering coefficient. In all simulations and computations, we uniformly set the network size N = 103 and γ = 1.0. Simulations begin with 20 randomly-chosen infectious seeds and the final infection size ρ is obtained by making average over 100 realizations of different initial infected nodes, which can reduce random fluctuation caused by the initial conditions. In order to compare the simulation and theoretical result, we compute values of max ( L ) and max ( L ) according to the power method. In the Fig. 1, we compare the monte-carlo simulations and theoretical results on the epidemic threshold, β c . For the stochastic simulation, when β > β c , ρ > 0 (active state); when β < β c , ρ = 0 (absorbing state). For the theoretical result, we check new parameters Rpqmf = max ( L ) + 1 for the pair QMF model (6) and Rqmf = max ( L ) + 1 for the QMF model (10). When Rpqmf or Rqmf = 1, the epidemic threshold is attained. Here, we take the response parameter α = 0.1 and α = 0.5. From Fig. 1, one can see that the simulation result is close to the theoretical prediction and the deviation between them becomes large when changing from the ER network to the BA network. This shows that the heterogeneity of network distribution can affect the accuracy of the pair quenched meanfield approach. As shown in [28], the prediction deviation will be reduced as the network size increases. Nevertheless, all the sub-
plots show that the pair QMF method is more accurate than the classic QMF method. Moreover, one can see that the epidemic threshold is enhanced when the response parameter α is raised, which indicates that the local behavioral response can inhibit the epidemic spreading for a given network. Besides, such influence on the epidemic threshold is not strong, especially for a BA scale-free network. This should be for the formulation of response function. Compared to the frequency-dependent function, the density-dependent response function has a small influence on the transmission rate [19]. 4. Conclusions and discussions We have considered how local behaviorial response affects the epidemic threshold in the SIS model on quenched networks. In order to analyze this problem, we adopted a more exact approach recently developed by Mata and Ferreira [28], i.e., the pair QMF approach, rather than the classical QMF approach. We were able to theoretically derive the mathematical formulation of the pair QMF model and show that the pair QMF approach holds under the standard moment closure. Furthermore, we built a mathematical model to analyze the effect of local behaviorial response on the epidemic threshold. The response parameter is included in the expression of max (L˜ ) and hence has an effect on the epidemic threshold. Let us take a random regular network (RRN) as an example to illustrate it. Suppose the network has a uniform degree ki = n for each node i. Its adjacency matrix has a largest eigenvalue n [40]. Therefore, the epidemic threshold for this case is given by −PQMF λ˜ RRN = c
1 1 1 1 + 2 α, n − 1 1 − αn n−1 n −n
showing that α has a significant influence on the epidemic threshold. In the present work, we assumed that the response function is linear and the epidemic information is density-dependent. We argue that this assumption is necessary for the model based on pair quenched mean-field approach. A general response function should
22
Q. Wu et al. / Chaos, Solitons and Fractals 96 (2017) 17–22
be more meaningful for the realistic case [17,19], so the theoretical approach used here needs further improved in the future . Acknowledgments This work was jointly supported by the NSFC grants (Nos: 61663015 and 61203153) and Natural Science Foundation of Jiangxi Province of China (20161BAB202051). Appendix A. The proof of pair approximation In order to close the equations, it is needed to use a pair approximation. Here, we would like to show that the standard pair approximation (2) holds for the random contact network with no clustering. In Eq. (2), Ai = {Xi (t ) = A}, B j = {X j (t ) = B} and Cl = {Xl (t ) = C } where A, B, C ∈ {0, 1}. Meanwhile, nodes i and l are two neighbors of node j. According to the definition of the conditional probability, Eq. (2) can be written as
P[Cl |Ai B j ]P[Ai B j ] =
P[Ai B j ]P[Cl |B j ]P[B j ] . P[B j ]
(A.1)
Note that when P[Ai B j ] = 0, the standard approximation (2) is true since Eq. (A.1) holds for this case. When P[Ai B j ] = 0, Eq. (A.1) is equivalent to
P[Cl |Ai B j ] = P[Cl |B j ], which means that the state of node l is uncorrelated to the state of its neighbor’s neighbor, e.g., node i. This condition holds approximately only for the static randomly connected network with no clustering. References [1] Pastor-Satorras R, Castellano C, Van Mieghem P, Vespignani A. Epidemic processes in complex networks. Rev Mod Phys 2015;87:925. [2] Nowzari C, Preciado VM, Pappas GJ. Analysis and control of epidemics: a survey of spreading processes on complex networks. IEEE Control Syst 2016;36(1):26–46. [3] Lindquist J, Ma J, van den Driessche P, Willeboordse FH. Effective degree network disease models. J Math Biol 2011;62(2):143–64. [4] Pastor-Satorras R, Vespignani A. Epidemic spreading in scale-free networks. Phys Rev Lett 20 01;86:320 0. [5] Mieghem PV, Omic J, Kooij RE. Virus spread in networks. IEEE ACM T Netw 2009;17:1. [6] Liu Y, Wei B, Du YX, Xiao FY, Deng Y. Identifying influential spreaders by weight degree centrality in complex networks. Chaos Solit Fract 2016;86:1–7. [7] Wang Z, Andrews MA, Wu ZX, Wang L, Bauch CT. Coupled disease-behavior dynamcis on complex networks: a review. Phys Life Rev 2015;15:1–29. [8] Shaw LB, Schwartz IB. Fluctuating epidemics on adaptive networks. Phys Rev E 2008;77(6):066101. [9] Granell C, Gomez S, Arenas A. Dynamical interplay between awareness and epidemic spreading in multiplex networks. Phys Rev Lett 2013;111(12): 128701. [10] Wang Z, Zhao DW, Wang L, Sun GQ, Jin Z. Immunity of multiplex networks via acquaintance vaccination. Europhys Lett 2015;112(4):48002. [11] Wang Z, Bauch CT, Bhattacharyya S. Statistical physics of vaccination. Phys Rep 2016;664:1–113.
[12] Sahneh FD, Chowdhury FN, Scoglio CM. On the existence of a threshold for preventive behavioral responses to suppress epidemic spreading. Sci Rep 2012;2(6):632. [13] Funk S, Gilad E, Jansen VA. Endemic disease, awareness, and local behavioural response. J Theor Biol 2010;264(2):501–9. [14] Zhang HF, Xie JR, Chen HS, Liu C, Small M. Impact of asymptomatic infection on coupled disease-behavior dynamics in complex networks. Europhys Lett 2016;114(3):38004. [15] Bagnoli F, Lió P, Sguanci L. Risk perception in epidemic modeling. Phys Rev E 2007;76(6):61904. [16] Wu QC, Fu XC, Small M, Xu XJ. The impact of awareness on epidemic spreading in network. Chaos 2012;22(1):013101. [17] Shang YL. Modeling epidemic spread with awareness and heterogeneous transmission rates in networks. J Biol Phys 2013;39(3):489–500. [18] Yan ZJ, Huang H, Chen YH, Pan YH. Identifying the direct risk source to contain epidemics more effectively. Phys Rev E 2016;93:012308. [19] Zhang HF, Xie JR, Tang M, Lal YC. Suppression of epidemic spreading in complex networks by local information based behavioral responses. Chaos 2014;24(4):043106. [20] Sun GQ, Zhang J, Song LP, Jin Z, Li BL. Pattern formation of a spatial predatorcprey system. Appl Math Comput 2012;218(22):11151–62. [21] Sun GQ, Wang SL, Ren Q, Jin Z, Wu YP. Effects of time delay and space on herbivore dynamics: linking inducible defenses of plants to herbivore outbreak. Sci Rep 2015;5:11246. [22] Sun GQ, Jusup M, Jin Z, Wang Y, Wang Z. Pattern transitions in spatial epidemics: mechanisms and emergent properties. Phys Life Rev 2016;19:43–73. [23] Ferreira SC, Castellano C, Pastor-Satorras R. Epidemic thresholds of the susceptible- infected-susceptible model on networks: a comparison of numerical and theoretical results. Phys Rev E 2012;86(4):041125. [24] Li C, van de Bovenkamp R, van Mieghem P. Susceptible-infected-susceptible model: a comparison of n-intertwined and heterogeneous mean-field approximations. Phys Rev E 2012;86(2):026116. [25] Wang Y, Chakrabarti D, Wang CX, et al. Epidemic spreading in real networks: an eigenvalue viewpoint. In: Proceedings of 22nd international symposium on reliable distributed systems. Pittsburgh, PA: Carnegie Mellon University; 2003. [26] Gómez S, Arenas A, Borge-Holthoefer J, Meloni S, Moreno Y. Discrete-time Markov chain approach to contact-based disease spreading in complex networks. Europhys Lett 2010;89:38009. [27] Cator E, van Mieghem P. Second-order mean-field susceptible-infected-susceptible epidemic threshold. Phys Rev E 2012;85(5):056111. [28] Mata AS, Ferreira SC. Pair quenched mean-field theory for the susceptible-infected-susceptible model on complex networks. Europhys Lett 2013;103:48003. [29] Wu QC, Zhang HF. Epidemic threshold of node-weighted susceptible-infected– susceptible models on networks. J Phys A 2016;49. 345601 (14pp) [30] Zhang XG, Sun GQ, Ma JL, Jin Z. Epidemic dynamics on semi-directed complex networks. Math Biosci 2013;246(2). 242-51 [31] Kostova TJ. Interplay of node connectivity and epidemic rates in the dynamics of epidemic networks. Differ Equ Appl 2009;15(4):415–28. [32] Mata AS, Ferreira RS, Ferreira SC. Heterogeneous pair-approximation for the contact process on complex networks. New J Phys 2014;16:053006. [33] Wu QC, Zhang HF, Zeng GH. Responsive immunization and intervention for infectious diseases in social networks. Chaos 2014;24:023108. [34] Luo XF, Zhang XG, Sun GQ, Jin Z. Epidemical dynamics of SIS pair approximation models on regular and radnom networks. Physica A 2014;410:144–53. [35] Eames KTD, Keeling MJ. Modeling dynamic and network heterogeneities in the spread of sexually transmitted diseases. Proc Natl Acad Sci U S A 2002;99(20):13330–5. [36] House T, Keeling MJ. Insights from unifying modern approximations to infections on networks. J R Soc Interface 2011;8(54):67–73. [37] Wu QC, Chen SF. Mean field theory of epidemic spreading with effective contacts on networks. Chaos Solit Fract 2015;81. 359-64 [38] Erdös P, Rényi A. On the evolution of random graphs. Publ Math Inst Hung Acad Sci 1960;5:17–61. [39] Barabási AL, Albert R. The emergence of scaling in random networks. Science 1999;286:509. [40] van Mieghem P. Graph spectra for complex networks. Cambridge: Cambridge University Press; 2011. p. 43–4.
Contents lists available at ScienceDirect
Chaos, Solitons and Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos
Epidemic spreading over quenched networks with local behavioral response Qingchu Wu a,∗, Shufang Chen b, Lingling Zha a a b
College of Mathematics and Information Science, Jiangxi Normal University, Nanchang, Jiangxi 330022, P. R. China College of Physics and Communication Electronics, Jiangxi Normal University, Nanchang, Jiangxi 330022, P. R. China
a r t i c l e
i n f o
Article history: Received 14 October 2016 Revised 20 December 2016 Accepted 8 January 2017
Keywords: Complex networks Epidemic spreading Pair approximation
a b s t r a c t We discuss the impact of local information-based behavioral response on epidemic spreading in social networks. By using a pair quenched mean-field approach developed by Mata and Ferreira [Europhys. Lett. 103 (2013) 48003], we derive a dynamical model governing the epidemic spreading over a random network with a linear response function and density-dependent epidemic information. A deterministic relation between the epidemic threshold and the response parameter is derived according to a quasi-static approximation method. It is found that local behavioral response will induce the extinction of the disease via rasing the epidemic threshold. Additionally, the theoretical result is supported by stochastic simulations on an Erdo¨ s–Rényi random network and a Baraba´ si–Albert scale-free network. Simulations show that the pair quenched mean-field approach is more accurate than the classical quenched mean-field approach. © 2017 Elsevier Ltd. All rights reserved.
1. Introduction The theory of complex network involves the network itself (structure and its dynamics), dynamics on networks and coupled interactions between network and dynamics evolving over the network [1]. As a typical process occurring on the network [2], disease transmission dynamics have drawn a wide attention of researchers from mathematical, computer, physical and biological communities [3–6]. Recently, many researches focus on the coupled disease-behavior dynamics [7]. Generally speaking, individual behavioral response means that individual adjusts its adaptive behavior by changing contact numbers or contact objects [8,9], vaccination decision-Erdo¨ s-Rényimaking [10,11] and dynamical parameters [12–14]. In reality, the epidemic information (denoted by x) can cause individuals to keep social distancing (by wearing protective masks, washing hands frequently, avoiding crowded public areas, or staying home from work or school or more creative precautions), which potentially results in the reduction of individual susceptibility (denoted by y) [13,16]. For a node with k contact number and s infected neighbors, if the epidemic information obtained by a node is determined by the infection fraction in its neighborhood, then
∗
Corresponding author: E-mail address: [email protected] (Q. Wu).
http://dx.doi.org/10.1016/j.chaos.2017.01.003 0960-0779/© 2017 Elsevier Ltd. All rights reserved.
x = s/k and we say it is density-dependent [15–18]; if the information is only related to the number of infected neighbors, then x = s and it is frequency-dependent [19]. To incorporate the information effect into the transmission model, it is needed to define a response function y = φ (x ), which satisfies φ (0 ) = 0, 0 ≤ φ (x ) ≤ 1, and φ (x) < 0. It is relevant to investigate the effect of individual behavioral response on the epidemic spreading under the network framework with the spatial diffusion [20–22]. In [16], we studied the susceptible-infected-susceptible (SIS) model with a response function of linear form, i.e., φ (x ) = 1 − α x(1 ≥ α ≥ 0 ). Meanwhile, the epidemic information is densitydependent. By using the heterogeneous mean-field (HMF) approach, we obtained the epidemic threshold above which the epidemic spreads over the whole network, depending on the response parameter α . Shang [17] analyzed the case of nonlinear form φ (x ) = 1 − α (s/k )m (m is a positive integer) and derived the epidemic threshold. Recently, Zhang et al. [19] proposed a response function of exponential form φ (x ) = (1 − α )x with the frequency-dependent epidemic information. Their theoretical analysis and simulations show that the epidemic threshold is strongly affected by individual response. So, the spread of an infectious disease over a network can be controlled by individual behaviorial response. Compared to the immunization and quarantine, individual response is a very economic strategy. All these researches are based on the HMF model, that has been shown to be exact only for the annealed network [23,24]. How to analyze this issue on the static or quenched network? [25] How to
18
Q. Wu et al. / Chaos, Solitons and Fractals 96 (2017) 17–22
study the interplay between the network structure and the behavioral response? For the SIS-type model, the microscopic Markov-chain (MMC) approach is frequently-used in the analysis of the epidemic spreading on the quenched network [26]. Van Mieghem et al. developed a continuous-time MMC approach and called it as the N-intertwined model [5], which is also referred to as the quenched mean-field (QMF) model [28]. Compared to the HMF modeling framework, the N-interwind model may be more adapted for the quenched network [24]. However, this model is not enough exact since it only applies one-order neighbor information. Therefore, a highorder mathematical formulation is required to provide a more exact threshold analysis [27]. Mata and Ferreira proposed a pair-type formulation of quenched mean-field model, referred to as the pair quenched mean-field model [28]. The pair QMF model makes use of the two-order neighbor information and gives a more exact prediction. More importantly, the epidemic threshold can be determined and analyzed by this model. In the present work, we would like to give a rigorous derivation for the pair QMF model, and show that the pair QMF model is adapted for the static randomly connected network with no clustering [28]. In addition, we analyze the impact of local behavioral response on the epidemic spreading over the static network by using the pair QMF model, which is still challenging for us since it involves the relation between the contact structure and local behavioral response.
Within this perspective, they built a system describing the SIS epidemic dynamics on the network G, which takes the form
2. An analysis framework
2.2. The derivation of the pair QMF equations
We assume that an infectious disease follows SIS dynamics on a given network with a size N, denoted by G = (V, E ). Generally, the contact network is weighted [29] or directed [30]. For simplicity, we assume that G is a unweighted and undirected network and hence it is completely determined by its adjacency matrix A = (ai j ): if node i links to node j in G, then ai j = 1; otherwise ai j = 0. In this model, each node may stay in either susceptible (S) or infected (I) state. During a infinitesimal time interval (t , t + t ], an infected node transmits the infectious disease to its susceptible neighbor with probability β t and meanwhile, it recovers and becomes susceptible again with probability γ t. As usual, we define λ = β /γ be the effective spreading rate [4,25]. The mathematical model for the static contact network can be built by the microscopic Markov-chain approach [26,31]. This approach focuses on the dynamic of the probability of each node i to be infected at time t, denoted by ρ i . Let Xi (t) ∈ {0, 1} denote the state of node i at time t [5]. If Xi (t ) = 0, node i is susceptible; if Xi (t ) = 1, node i is infectious. Then ρi = P[Xi (t ) = 1] and 1 − ρi = P[Xi (t ) = 0], where P[A] represents the probability of event A occurring.
We will apply the total probability formula to derive the pair QMF equations for the SIS epidemic model (1). For the sake of the following analysis, we first present a lemma.
2.1. The pair quenched mean-field model Mata and Ferreira [28] proposed a pair-type formulation of quenched mean-field model, i.e., the pair QMF model. In order to develop the pair QMF theory, they introduced a series of simple notations: [Ai ] = P[Xi (t ) = A], which is the probability that node i is in the state A; [Ai , B j ] = P[Xi (t ) = A, X j (t ) = B], which is the probability that nodes i and j are in states A and B, respectively; [Ai , Bj , Cl ] is the extension to three nodes. In addition, they also defined specific pair-type variables: ωi j = [0i , 0 j ], φi j = [0i , 1 j ], φ¯ i j = [1i , 0 j ], ψi j = [1i , 1 j ]. It is easy to know that pair-type variables must satisfy
ψi j = ρ j − φi j , φ¯ i j = ρi − ψi j = ρi − ρ j + φi j , ωi j = 1 − ρi − φi j .
d ρi = −γ ρi + β φi j a i j dt N
j=1
d φi j = −γ φi j − βφi j + γ ψi j dt +β [0i , 0 j , 1l ] − β [1l , 0i , 1 j ] l∈N ( j ) l=i
(1)
l∈N (i ) l= j
where N (i ) denotes the neighborhood of node i. In order to close the above equations, the authors used the standard approximation [28,32]
[Ai , B j , Cl ]
[Ai , B j ][B j , Cl ] . [B j ]
(2)
In Appendix A, we prove that this approximation holds for the random network with no clustering. After performing a quasi-static approximation for t → ∞, the authors found that the epidemic threshold can be determined by the largest eigenvalue of the Jacobian matrix L, which is given by
Li j = −
γ+
β 2 ki β ( 2γ + β ) ai j . δi j + 2β + 2γ 2β + 2γ
(3)
Here, δ ij is the Kronecker symbol and ki is the degree of node i.
Lemma 1. For each a subset of nodes in G, V∗ (i.e., V∗ ⊂ V), we assume that a probability of each node i ∈ V∗ to be infected, denoted by ρ i is given, then the number of infected nodes in this subset is a stochastic variable ξ ∈ [0, |V∗ |] (|V∗ | denotes the number of the elements in set V∗ ) and its first and second order moment satisfy
E[ξ ] =
ρi , E [ ξ 2 ] =
i∈V ∗
ρi +
i∈V ∗
ρi 1 ρi 2 ,
i1 =i2
respectively. This lemma is easily proved by the mathematical induction. In fact, the first order moment formula has been proved in [31,33]. The second order formula can be justified by a simple example when |V ∗ | = 3. At this time, E[ξ 2 ] = 12 × [ρ1 (1 − ρ2 )(1 − ρ3 ) + ρ2 (1−ρ1 )(1 − 3 ) + ρ3 (1−ρ1 )(1−ρ2 )] + 22 × [ρ1 ρ2 (1 − ρ3 ) + ρ1 ρ3 (1 − ρ2 ) + ρ2 ρ3 (1 − ρ1 )] + 32 × ρ1 ρ2 ρ3 = ρ1 + ρ2 + ρ3 + 2(ρ1 ρ2 + ρ2 ρ3 + ρ1 ρ3 ), which accords with the conclusion from Lemma 1. In the following, we derive system (1). By the total probability formula, we have
P[Xi (t + t ) = 1] = P[Xi (t ) = 1]P[Xi (t + t ) = 1|Xi (t ) = 1] + P[Xi (t ) = 0]P[Xi (t + t ) = 1|Xi (t ) = 0].
(4)
Denote kinf (i) be the number of infected nodes in the neighborhood of node i. Note that
P[Xi (t ) = 0]P[Xi (t + t ) = 1|Xi (t ) = 0] = P[Xi (t + t ) = 1, Xi (t ) = 0] =
ki
P[Xi (t + t ) = 1, Xi (t ) = 0, kinf (i ) = s]
s=0
=
ki s=0
P[Xi (t + t ) = 1|Xi (t ) = 0, kinf (i ) = s]
Q. Wu et al. / Chaos, Solitons and Fractals 96 (2017) 17–22
P[Xi (t + t ) = 0, X j (t + t ) = 1|Xi (t ) = 1, X j (t ) = 0] = o(t ).
P[Xi (t ) = 0, kinf (i ) = s] ki
=
[1 − (1 − β t )s ]P[Xi (t ) = 0, kinf (i ) = s]
s=0
= E pa [1 − (1 − β t )s ] Then, we have
P[Xi (t + t ) = 1] = P[Xi (t ) = 1]P[Xi (t + t ) = 1|Xi (t ) = 1] +
ki
[1 − (1 − β t )s ]P[Xi (t ) = 0, kinf (i ) = s]
s=0
= P[Xi (t ) = 1](1 − γ ) +
ki
[1 − (1 − β t )s ]P[Xi (t ) = 0, kinf (i ) = s]
s=0
Thus, subtracting P[Xi (t ) = 1] and dividing by t and letting t → 0, we get
d ρi = −γ ρi + β dt
ki
sP[Xi (t ) = 0, kinf (i ) = s].
s=0
Xi (t ) = 0, X j (t ) = 0, kinf ( j ) = s]
β t P[Xi (t ) = 0, X j (t ) = 1, kinf (i ) − 1 = s]
ki −1
(5)
P[Xi (t + t ) = 1, X j (t + t ) = 1, Xi (t ) = 0, X j (t ) = 1,
s=0
kinf (i ) = s] ki −1
=
P[Xi (t + t ) = 1, X j (t + t ) = 1|Xi (t ) = 0, X j (t ) = 1,
s=0
kinf (i ) = s] P[Xi (t ) = 0, X j (t ) = 1, kinf (i ) = s] ki −1
=
[1 − (1 − β t )s ]P[Xi (t ) = 0, X j (t ) = 1, kinf (i ) = s]
s=0
k j −1
P[Xi (t + t ) = 0, X j (t + t ) = 1|Xi (t ) = 0,
s=0
X j (t ) = 0, kinf ( j ) = s] P[Xi (t ) = 0, X j (t ) = 0, kinf ( j ) = s] k j −1
[1 − (1 − β t )s ]P[Xi (t ) = 0, X j (t ) = 0, kinf ( j ) = s]
s=0
s = Etr 1 [1 − (1 − β t ) ],
and
P[Xi (t + t ) = 0, X j (t + t ) = 1|Xi (t ) = 1, X j (t ) = 1] = γ t + o(t ), and
ki −1
=
It is needed to give the specific expression of the four terms on the right-hand-side (RHS) of Eq. (5). For the former three terms, we have
P[Xi (t + t ) = 0, X j (t + t ) = 1,
P[Xi (t ) = 0, X j (t ) = 1, kinf (i ) = s]
P[Xi (t + t ) = 1, X j (t + t ) = 1, Xi (t ) = 0, X j (t ) = 1]
+ P[Xi (t + t ) = 0, X j (t + t ) = 1|Xi (t ) = 0, X j (t ) = 1]
k j −1
kinf (i ) − 1 = s]
For part (ii), we turn to compute the probability P[Xi (t + t ) = 1, X j (t + t ) = 1, Xi (t ) = 0, X j (t ) = 1], which obeys
P[Xi (t ) = 1, X j (t ) = 1]
P[Xi (t + t ) = 0, X j (t + t ) = 1, Xi (t ) = 0, X j (t ) = 0]
P[Xi (t + t ) = 1, X j (t + t ) = 1|Xi (t ) = 0, X j (t ) = 1,
s=0
+ P[Xi (t + t ) = 0, X j (t + t ) = 1|Xi (t ) = 1, X j (t ) = 1]
P[Xi (t ) = 0, X j (t ) = 1]
= Etr 2 [β t] = β t P[Xi (t ) = 0, X j (t ) = 1].
+ P[Xi (t + t ) = 0, X j (t + t ) = 1|Xi (t ) = 1, X j (t ) = 0]
P[Xi (t + t ) = 1, X j (t + t ) = 1, Xi (t ) = 0, X j (t ) = 1,
s=0
P[Xi (t ) = 1, X j (t ) = 0]
=
ki −1
=
=
= P[Xi (t + t ) = 0, X j (t + t ) = 1, Xi (t ) = 0, X j (t ) = 0]
ki −1
kinf (i ) − 1 = s]
s=0
P[Xi (t + t ) = 0, X j (t + t ) = 1]
=
P[Xi (t + t ) = 1, X j (t + t ) = 1, Xi (t ) = 0, X j (t ) = 1]
s=0
Following Lemma 1, s=0 sP[Xi (t ) = 0, kinf (i ) = s] = j∈N (i ) P[Xi (t ) = 0, X j (t ) = 1]. So the first equation of Eq. (1) is derived. From this equation, one finds that the term φ ij has no influence on the dynamic of ρ i when node i does not connect to node j. Hence, it is required for us to further analyze the dynamic of φ ij as a pair of connected node (i, j). By the total probability formula, we have
While, the computation of the fourth term in the RHS of Eq. (5) is not straightforward. We notice that P[Xi (t + t ) = 0, X j (t + t ) = 1|Xi (t ) = 0, X j (t ) = 1] = 1−P[Xi (t + t ) = 0, X j (t + t ) = 0|Xi (t ) = 0, X j (t ) = 1] − P[Xi (t + t ) = 1, X j (t + t ) = 0|Xi (t ) = 0, X j (t ) = 1] − P[Xi (t + t ) = 1, X j (t + t ) = 1|Xi (t ) = 0, X j (t )=1]. Clearly, P[Xi (t + t ) = 0, X j (t + t ) = 0|Xi (t ) = 0, X j (t ) = 1] = γ t and P[Xi (t + t ) = 1, X j (t + t ) = 0|Xi (t ) = 0, X j (t ) = 1] = o(t ). It is a little complicated for the computation of P[Xi (t + t ) = 1, X j (t + t ) = 1|Xi (t ) = 0, X j (t ) = 1] because we do not know about the number of infected neighbors of node i. To solve this, we decompose it to two parts: (i) node i is infected by node j since j is infectious; (ii) node i is infected by node l (l = j. For part (i), it is easy to get that P[Xi (t + t ) = 1, X j (t + t ) = 1|Xi (t ) = 0, X j (t ) = 1] = β t. However, this also can be obtained by using the property of the joint probability:
=
k i
=
19
s = Etr 3 [1 − (1 − β t ) ].
Finally, subtracting P[Xi (t ) = 0, X j (t ) = 1] and dividing by t and letting t → 0 in Eq. (5) and using Lemma 1, one can obtain the second equation of Eq. (1). It is worthy noting that the derivation of Eq. (1) does not need any external assumption, so it is true for any network structure. 3. The local behavioral response model 3.1. The model The objective of the present work is to investigate the effect of the local epidemic information on the epidemic dynamic of SIS model. In our model, it is assumed that the response function is a linear formulation with respect to the epidemic information and is also density-dependent [16]. It suffices to replace β
20
Q. Wu et al. / Chaos, Solitons and Fractals 96 (2017) 17–22
by β [1 − α kinf (i )/ki ] for susceptible node i. We apply the analysis framework as mentioned in Section 2 to derive corresponding pair QMF equations. Note that
E [1 − (1 − β (1 − α s/ki )t ) ] t pa
lim
t→0
s
βα
= β E pa [s] −
E pa [s2 ] ki α βα =β 1− [0i , 1 j ] − ki ki j∈N (i )
N α
=β 1− lim
t→0
ki
φi j a i j − β
j=1
ki
=β
=β 1−
−
α
βα
[0i , 1 j1 ][0i , 1 j2 ]
d α ρi = −γ ρi + β 1 − dt ki
j1 , j2 ∈N (i ), j1 = j2
φi j 1 φi j 2 a i j 1 a i j 2 ,
kj
+β 1 −
[0i , 0 j , 1l ] −
l∈N ( j ) l=i
βα kj
[0i , 0 j , 1l1 ][0i , 0 j , 1l2 ],
l1 ,l2 ∈N ( j ) l1 ,l2 =i,l1 =l2
Etr βα tr 2 [β (1 − α (s + 1 )/ki )t] = β Etr E [s] 2 [1] − t→0 t ki 2 α βα =β 1− [1 , 0i , 1 j ] φi j − ki ki l∈N (i) l lim
l= j
and
lim
t→0
s Etr 3 [1 − (1 − β (1 − α s/ki )t ) ] t
= β Etr 3 [s] −
Etr [s2 ] ki 3 α βα =β 1− [1 , 0i , 1 j ] − ki l∈N (i) l ki l= j
d φi j = dt
β 1− −
βα kj
α kj
N
[1l1 , 0i , 1 j ][1l2 , 0i , 1 j ].
l1 ,l2 ∈N (i ) l1 ,l2 = j,l1 =l2
φi j a i j − β
j=1
+
βα ki
ki
φi j 1 φi j 2 a i j 1 a i j 2
ki
j
j
2γ + β (1 − kα ) + β (1 − kα )
.
(8)
j
ui j =
β (1 − α /k j ) , 2γ + β (1 − α /ki ) + β (1 − α /k j )
vi j =
2γ + β (1 − α /k j ) . 2γ + β (1 − α /ki ) + β (1 − α /k j )
After plugging (8) in (7), we obtain
d ρi = L˜i j ρ j , dt
j1 = j2
where the system Jacobian is given by
⎛ βα ⎝ ki
1+
[1l , 0i , 1 j ]⎠
l∈N (i ) l= j
[1l , 0i , 1 j ]
l∈N (i ) l= j
[1l1 , 0i , 1 j ][1l2 , 0i , 1 j ]
l1 ,l2 ∈N (i ) l1 ,l2 = j,l1 =l2
ωi j φ jl φ¯ li φi j , [1l , 0i , 1 j ] = . 1 − ρj 1 − ρi
γ +β 1−
N α
ki
ai j ui j
α δi j + β 1 − ai j vi j . ki
j=1
˜ PQMF , is deterAs a result, the epidemic threshold, denoted by λ c mined by the equation max ( L ) = 0. Here, max (M) represents the largest eigenvalue of the matrix M. So, the outbreak condition can be seen to depend only of the values of α , A and β /γ . Finally, we analyzed the local response model by taking an one-vertex approximation, i.e., P[Xi (t ) = 0, X j (t ) = 1] = P[Xi (t ) = 0]P[X j (t ) = 1], or φi j (1 − ρi )ρ j . At this time, the epidemic spreading is simply governed by
⎞
(9)
j=1
[0i , 0 j , 1l1 ][0i , 0 j , 1l2 ]
α
[2γ + β (1 − kα )]ρ j − β (1 − kα )ρi
Li j = −
where
[0i , 0 j , 1l ] =
α
[0i , 0 j , 1l ]
+ γ ψi j − γ φi j − βφi j +
(7)
N
l1 ,l2 ∈N ( j ) l1 ,l2 =i,l1 =l2
−β 1 −
φ jl (a jl − δil )
l=1
For clarity, let
l∈N ( j ) l=i
kj
N
φi j + γ ρ j
j
Upon substituting these into the Eqs. (4) and (5), we write down equations for the changes in the SIS model with local behavioral response
d α ρi = −γ ρi + β 1 − dt ki
α
φi j a i j
j=1
However, it seems hard to derive the epidemic threshold for system (7) according to the stability of its Jacobin matrix [34], due to its complicated formulation. Indeed, as for the traditional pairapproximation model on heterogeneous networks [35,36], it is not easy to explicitly derive analytical expressions of the outbreak condition. However, following [28,32,37], a quasi-static approximation can be applied into our model to solve this problem. For t → ∞, we let dρ i /dt ࣃ 0 and dφ i, j /dt ࣃ 0 when the system is around the zero solution. Hence, we have β (1 − kα ) j N N l=1 φ jl a jl = γ ρ j and l=1 φ jl δil = φ ji . Then, substituting these equations in Eq. (7) leads to
φi j =
βα
N
d α φi j = − 2 γ + β 1 − dt ki
j1 = j2
2 Etr 1 [s ]
kj
j
ψi j + φi j = ρ j , we have
s Etr 1 [1 − (1 − β (1 − α s/k j )t ) ] t
Etr 1 [s]
In the epidemiology, it is critical to determine the epidemic threshold, by which one can judge whether the infectious disease spreads in a population. To this end, we write a linear form of the model around the fixed point ρi = φi j = ψi j = 0. At this time, N ωi j φ jl N l∈N ( j ) [0i , 0 j , 1l ] = l=1 1−ρ (a jl − δil ) l=1 φ jl (a jl − δil ). Since l=i
α
3.2. Epidemic threshold
(6)
d ρi = −γ ρi dt + β ( 1 − ρi )
1−
N α
ki
j=1
ρ j ai j −
α ki
ρ j1 ρ j2 ai j1 ai j2 .
j1 = j2
(10) Denote Li j be the Jacobian matrix of Eq. (10). When its largest ˆ QMF ) eigenvalue is null, the epidemic threshold (denoted by λ c
Q. Wu et al. / Chaos, Solitons and Fractals 96 (2017) 17–22
21
L + 1 and final infection size ρ as functions of the infection rate β . Rpqmf = 1 and Rqmf = 1 imply that max (L˜ ) = 0 Fig. 1. Illustrations of Rpqmf = max ( L ) + 1, Rqmf = max L = 0, respectively. They represent the critical levels of epidemic outbreak. and max
satisfies
λˆ QMF = c
1
max (H )
,
(11)
where H = [hi j ]N×N , hi j = ai j − α ai j /ki . 3.3. Simulations To test above argument, we perform stochastic simulations over both two typical contact networks. One is a Baraba´ si–Albert (BA) scale-free network with the average degree k = 6 [39], which has a highly heterogeneous degree distribution (power-law distribution) and a small clustering coefficient. The other is an Erdo¨ s–Rényi (ER) random network with connecting probability p = 0.006 [38], which has a homogeneous degree distribution (Poisson distribution) and a small clustering coefficient. In all simulations and computations, we uniformly set the network size N = 103 and γ = 1.0. Simulations begin with 20 randomly-chosen infectious seeds and the final infection size ρ is obtained by making average over 100 realizations of different initial infected nodes, which can reduce random fluctuation caused by the initial conditions. In order to compare the simulation and theoretical result, we compute values of max ( L ) and max ( L ) according to the power method. In the Fig. 1, we compare the monte-carlo simulations and theoretical results on the epidemic threshold, β c . For the stochastic simulation, when β > β c , ρ > 0 (active state); when β < β c , ρ = 0 (absorbing state). For the theoretical result, we check new parameters Rpqmf = max ( L ) + 1 for the pair QMF model (6) and Rqmf = max ( L ) + 1 for the QMF model (10). When Rpqmf or Rqmf = 1, the epidemic threshold is attained. Here, we take the response parameter α = 0.1 and α = 0.5. From Fig. 1, one can see that the simulation result is close to the theoretical prediction and the deviation between them becomes large when changing from the ER network to the BA network. This shows that the heterogeneity of network distribution can affect the accuracy of the pair quenched meanfield approach. As shown in [28], the prediction deviation will be reduced as the network size increases. Nevertheless, all the sub-
plots show that the pair QMF method is more accurate than the classic QMF method. Moreover, one can see that the epidemic threshold is enhanced when the response parameter α is raised, which indicates that the local behavioral response can inhibit the epidemic spreading for a given network. Besides, such influence on the epidemic threshold is not strong, especially for a BA scale-free network. This should be for the formulation of response function. Compared to the frequency-dependent function, the density-dependent response function has a small influence on the transmission rate [19]. 4. Conclusions and discussions We have considered how local behaviorial response affects the epidemic threshold in the SIS model on quenched networks. In order to analyze this problem, we adopted a more exact approach recently developed by Mata and Ferreira [28], i.e., the pair QMF approach, rather than the classical QMF approach. We were able to theoretically derive the mathematical formulation of the pair QMF model and show that the pair QMF approach holds under the standard moment closure. Furthermore, we built a mathematical model to analyze the effect of local behaviorial response on the epidemic threshold. The response parameter is included in the expression of max (L˜ ) and hence has an effect on the epidemic threshold. Let us take a random regular network (RRN) as an example to illustrate it. Suppose the network has a uniform degree ki = n for each node i. Its adjacency matrix has a largest eigenvalue n [40]. Therefore, the epidemic threshold for this case is given by −PQMF λ˜ RRN = c
1 1 1 1 + 2 α, n − 1 1 − αn n−1 n −n
showing that α has a significant influence on the epidemic threshold. In the present work, we assumed that the response function is linear and the epidemic information is density-dependent. We argue that this assumption is necessary for the model based on pair quenched mean-field approach. A general response function should
22
Q. Wu et al. / Chaos, Solitons and Fractals 96 (2017) 17–22
be more meaningful for the realistic case [17,19], so the theoretical approach used here needs further improved in the future . Acknowledgments This work was jointly supported by the NSFC grants (Nos: 61663015 and 61203153) and Natural Science Foundation of Jiangxi Province of China (20161BAB202051). Appendix A. The proof of pair approximation In order to close the equations, it is needed to use a pair approximation. Here, we would like to show that the standard pair approximation (2) holds for the random contact network with no clustering. In Eq. (2), Ai = {Xi (t ) = A}, B j = {X j (t ) = B} and Cl = {Xl (t ) = C } where A, B, C ∈ {0, 1}. Meanwhile, nodes i and l are two neighbors of node j. According to the definition of the conditional probability, Eq. (2) can be written as
P[Cl |Ai B j ]P[Ai B j ] =
P[Ai B j ]P[Cl |B j ]P[B j ] . P[B j ]
(A.1)
Note that when P[Ai B j ] = 0, the standard approximation (2) is true since Eq. (A.1) holds for this case. When P[Ai B j ] = 0, Eq. (A.1) is equivalent to
P[Cl |Ai B j ] = P[Cl |B j ], which means that the state of node l is uncorrelated to the state of its neighbor’s neighbor, e.g., node i. This condition holds approximately only for the static randomly connected network with no clustering. References [1] Pastor-Satorras R, Castellano C, Van Mieghem P, Vespignani A. Epidemic processes in complex networks. Rev Mod Phys 2015;87:925. [2] Nowzari C, Preciado VM, Pappas GJ. Analysis and control of epidemics: a survey of spreading processes on complex networks. IEEE Control Syst 2016;36(1):26–46. [3] Lindquist J, Ma J, van den Driessche P, Willeboordse FH. Effective degree network disease models. J Math Biol 2011;62(2):143–64. [4] Pastor-Satorras R, Vespignani A. Epidemic spreading in scale-free networks. Phys Rev Lett 20 01;86:320 0. [5] Mieghem PV, Omic J, Kooij RE. Virus spread in networks. IEEE ACM T Netw 2009;17:1. [6] Liu Y, Wei B, Du YX, Xiao FY, Deng Y. Identifying influential spreaders by weight degree centrality in complex networks. Chaos Solit Fract 2016;86:1–7. [7] Wang Z, Andrews MA, Wu ZX, Wang L, Bauch CT. Coupled disease-behavior dynamcis on complex networks: a review. Phys Life Rev 2015;15:1–29. [8] Shaw LB, Schwartz IB. Fluctuating epidemics on adaptive networks. Phys Rev E 2008;77(6):066101. [9] Granell C, Gomez S, Arenas A. Dynamical interplay between awareness and epidemic spreading in multiplex networks. Phys Rev Lett 2013;111(12): 128701. [10] Wang Z, Zhao DW, Wang L, Sun GQ, Jin Z. Immunity of multiplex networks via acquaintance vaccination. Europhys Lett 2015;112(4):48002. [11] Wang Z, Bauch CT, Bhattacharyya S. Statistical physics of vaccination. Phys Rep 2016;664:1–113.
[12] Sahneh FD, Chowdhury FN, Scoglio CM. On the existence of a threshold for preventive behavioral responses to suppress epidemic spreading. Sci Rep 2012;2(6):632. [13] Funk S, Gilad E, Jansen VA. Endemic disease, awareness, and local behavioural response. J Theor Biol 2010;264(2):501–9. [14] Zhang HF, Xie JR, Chen HS, Liu C, Small M. Impact of asymptomatic infection on coupled disease-behavior dynamics in complex networks. Europhys Lett 2016;114(3):38004. [15] Bagnoli F, Lió P, Sguanci L. Risk perception in epidemic modeling. Phys Rev E 2007;76(6):61904. [16] Wu QC, Fu XC, Small M, Xu XJ. The impact of awareness on epidemic spreading in network. Chaos 2012;22(1):013101. [17] Shang YL. Modeling epidemic spread with awareness and heterogeneous transmission rates in networks. J Biol Phys 2013;39(3):489–500. [18] Yan ZJ, Huang H, Chen YH, Pan YH. Identifying the direct risk source to contain epidemics more effectively. Phys Rev E 2016;93:012308. [19] Zhang HF, Xie JR, Tang M, Lal YC. Suppression of epidemic spreading in complex networks by local information based behavioral responses. Chaos 2014;24(4):043106. [20] Sun GQ, Zhang J, Song LP, Jin Z, Li BL. Pattern formation of a spatial predatorcprey system. Appl Math Comput 2012;218(22):11151–62. [21] Sun GQ, Wang SL, Ren Q, Jin Z, Wu YP. Effects of time delay and space on herbivore dynamics: linking inducible defenses of plants to herbivore outbreak. Sci Rep 2015;5:11246. [22] Sun GQ, Jusup M, Jin Z, Wang Y, Wang Z. Pattern transitions in spatial epidemics: mechanisms and emergent properties. Phys Life Rev 2016;19:43–73. [23] Ferreira SC, Castellano C, Pastor-Satorras R. Epidemic thresholds of the susceptible- infected-susceptible model on networks: a comparison of numerical and theoretical results. Phys Rev E 2012;86(4):041125. [24] Li C, van de Bovenkamp R, van Mieghem P. Susceptible-infected-susceptible model: a comparison of n-intertwined and heterogeneous mean-field approximations. Phys Rev E 2012;86(2):026116. [25] Wang Y, Chakrabarti D, Wang CX, et al. Epidemic spreading in real networks: an eigenvalue viewpoint. In: Proceedings of 22nd international symposium on reliable distributed systems. Pittsburgh, PA: Carnegie Mellon University; 2003. [26] Gómez S, Arenas A, Borge-Holthoefer J, Meloni S, Moreno Y. Discrete-time Markov chain approach to contact-based disease spreading in complex networks. Europhys Lett 2010;89:38009. [27] Cator E, van Mieghem P. Second-order mean-field susceptible-infected-susceptible epidemic threshold. Phys Rev E 2012;85(5):056111. [28] Mata AS, Ferreira SC. Pair quenched mean-field theory for the susceptible-infected-susceptible model on complex networks. Europhys Lett 2013;103:48003. [29] Wu QC, Zhang HF. Epidemic threshold of node-weighted susceptible-infected– susceptible models on networks. J Phys A 2016;49. 345601 (14pp) [30] Zhang XG, Sun GQ, Ma JL, Jin Z. Epidemic dynamics on semi-directed complex networks. Math Biosci 2013;246(2). 242-51 [31] Kostova TJ. Interplay of node connectivity and epidemic rates in the dynamics of epidemic networks. Differ Equ Appl 2009;15(4):415–28. [32] Mata AS, Ferreira RS, Ferreira SC. Heterogeneous pair-approximation for the contact process on complex networks. New J Phys 2014;16:053006. [33] Wu QC, Zhang HF, Zeng GH. Responsive immunization and intervention for infectious diseases in social networks. Chaos 2014;24:023108. [34] Luo XF, Zhang XG, Sun GQ, Jin Z. Epidemical dynamics of SIS pair approximation models on regular and radnom networks. Physica A 2014;410:144–53. [35] Eames KTD, Keeling MJ. Modeling dynamic and network heterogeneities in the spread of sexually transmitted diseases. Proc Natl Acad Sci U S A 2002;99(20):13330–5. [36] House T, Keeling MJ. Insights from unifying modern approximations to infections on networks. J R Soc Interface 2011;8(54):67–73. [37] Wu QC, Chen SF. Mean field theory of epidemic spreading with effective contacts on networks. Chaos Solit Fract 2015;81. 359-64 [38] Erdös P, Rényi A. On the evolution of random graphs. Publ Math Inst Hung Acad Sci 1960;5:17–61. [39] Barabási AL, Albert R. The emergence of scaling in random networks. Science 1999;286:509. [40] van Mieghem P. Graph spectra for complex networks. Cambridge: Cambridge University Press; 2011. p. 43–4.