The spreading of an infectious disease can trigger human behavior responses to the disease, which in turn plays a crucial role on the spreading of epidemic. In this study, to illustrate the impacts of the human behavioral responses, a new class of individuals, SF, is introduced to the classical susceptible-infected-recovered model. In the model, SF state represents that susceptible individuals who take self-initiate protective measures to lower the probability of being infected, and a susceptible individual may go to SF state with a response rate when contacting an infectious neighbor. Via the percolation method, the theoretical formulas for the epidemic threshold as well as the prevalence of epidemic are derived. Our finding indicates that, with the increasing of the response rate, the epidemic threshold is enhanced and the prevalence of epidemic is reduced. The analytical results are also verified by the numerical simulations. In addition, we demonstrate that, because the mean field method neglects the dynamic correlations, a wrong result based on the mean field method is obtained—the epidemic threshold is not related to the response rate, i.e., the additional SF state has no impact on the epidemic threshold.

Human behavioral responses to the spreading of epidemics have been recognized to have great influence on the epidemic dynamics. Therefore, it is very important to incorporate human behaviors into epidemiological models, which could improve models' utility in reflecting the reality and studying corresponding controlling measures. However, analytically grounded approaches to the problem of interacting effect between epidemic dynamics and human behavioral responses are still lacking so far, and what can help us better understand the impacts of human behavioral responses. In this work, by introducing a new SF state into classical susceptible-infected-recovered (SIR) model to mimic the situation that, when susceptible individuals are aware of the risk of infection, they may take self-protective measures to lower the probability of being infected. We then derive theoretical formulas based on the percolation method for the epidemic threshold and the prevalence of epidemic. Our results indicate that the introduction of the new SF state can significantly enhance the epidemic threshold and reduce the prevalence of epidemic. It is worth mentioning that a wrong result may be obtained when using the traditional mean field method—the epidemic threshold is not altered by the additional SF state. The result highlights that if the effects of human behavioral responses are ignored in mathematical modeling, the obtained results cannot really reflect the spreading mechanism of epidemics among human population.

Since network models can well describe the spreading of infectious disease among populations, many efforts have been devoted to studying this field.1,2 At first, researchers mainly paid attention to analyze the impact of the network structure on epidemic spreading and the control strategies, for example, how the network topology affects the epidemic threshold and the prevalence of epidemic,3–8 or how to design an effective immunization strategy to control the outbreaks of epidemics.9,10 In reality, outbreaks of epidemics can trigger spontaneous behavioral responses of individuals to take preventive measures, which in turn alters the epidemic dynamics and affects the disease transmission process.11–17 Thus, recently, some studies have made attempts to evaluate the impact of the human behaviors on epidemic dynamics. For instance, Funk et al.18 studied the impacts of awareness spread on both epidemic threshold and prevalence, and they found that, in a well-mixed population, spread of awareness can reduce the prevalence of epidemic but does not tend to affect the epidemic threshold, yet the epidemic threshold is altered when considering on social networks; Sahneh et al. considered a Susceptible-Alter-Infected-Susceptible (SAIS) model,19 and they found that the way of behavioral response can enhance the epidemic threshold; Meloni et al. constructed a meta-population model incorporating several scenarios of self-initiated behavioral changes into the mobility patterns of individuals, and they found that such behavioral changes do not alter the epidemic threshold, but may increase rather than decrease the prevalence of epidemic.20 Meanwhile, in Ref. 21, by designing the transmission rate of epidemic as a function of the local infected density or the global infected density, Wu et al. investigated the effect of such a behavioral response on the epidemic threshold.

One fact is that the infectious neighbors can infect a susceptible individual, also triggering the awareness of the susceptible individual.19,22 In view of this, in Ref. 23, Perra et al. introduced a new class of individuals, SF, that represents susceptible people who self-initiate behavioral changes that lead to a reduction in the transmissibility of the infectious disease, into the classical SIR model, and they found that such a model (SSFIR) can induce a rich phase space with multiple epidemic peaks and tipping points. However, the network structure was not incorporated into these models. As we know, when the model is considered within the network based framework, new theoretical tools should be used and new phenomena may be observed. In view of this, we incorporate the network structure into the SSFIR model23 to investigate its spreading dynamics. In the model, when contacting an infectious neighbor, susceptible individuals may be infected (from S state I state) with a transmission rate, or a behavioral response may be triggered (from S state to SF state) with a response rate. We provide a theoretical formula for the epidemic threshold as well as the prevalence of epidemic via the percolation method,5 our results show that the introduction of SF class can enhance the epidemic threshold and reduce the prevalence of epidemic. We also demonstrate that a wrong result can be obtained—the introduction of SF class cannot alter the epidemic threshold when using mean field method to such a model. The reasons are presented in Sec. V.

For the classical SIR model on complex network, where each node on network can be in one of three states: susceptible (S), infected (I), or recovered (R). The transmission rate along each SI link is β, and an infected node can enter R state with a recovery rate μ. To reflect the fact that, upon observation of infection, susceptible individuals may adopt protective measures to lower their infection risk, a new class, denoted by SF, is introduced into the original SIR model, we use SSFIR model to denote the modified SIR model in this study. In the model, when an S node contacts an I neighbor, besides the probability of being infected, the S node can go to SF state with a response rate βF0. The transmission rate for the SF nodes is lower; thus, we assume the transmission rate for SF nodes is γβ, where 0γ<1 is a discount factor.

The SSFIR model is described by the four following reactions and the associated rates:

S+Iβ2I,
(1)
S+IβFSF+I,
(2)
SF+Iγβ2I,
(3)
IμR.
(4)

Note that the SSFIR model returns to the SIR model once βF=0, and the SF state corresponds to fully vaccinated state when γ = 0.

In our model, during a sufficiently small time interval Δt, the transition rates of an SI edge becoming an II, SFI, and SR edge are β, βF, and μ, respectively. As a result, the probabilities of an SI edge becoming an II and SFI edge are given as T1=ββ+βF+μ and T2=βFβ+βF+μ, respectively. Similarly, since the transition rate of an SFI edge becoming an II and SFR during a sufficiently small time interval Δt are γβ and μ, the probability of an SFI edge becoming an II edge is T3=γβγβ+μ.24 

To analyze our proposed model, we first define “externally infected neighbor” (EIN) for any node. For node i, a neighbor j is an EIN, meaning that j is infected by its neighbors other than i (i.e., j is infected even though node i is deleted from the networks), which is the basic assumption of the cavity theory in statistical physics. Note that this method is suitable for the networks with negligible number of loops as the network size is sufficiently large25). The probability of neighbor j being an EIN of i is defined as θ; then, for the a node with degree k, the probability of having m EINs is given as

p(m|k)=(mk)θm(1θ)km.
(5)

Let p(R|m) be the probability of i being infected when it has number m of EIN. To calculate such a probability, we need to recognize that, in our model, an S node can become an I node through two ways: (a) the S node is directly infected; or (b) the S node first goes to SF state and then is infected. To facilitate the analysis, we approximately assume that the impacts of i's infected neighbors on node i happen in a non-overlapping order, i.e., they play roles on node i one by one.

For case (a), the probability of node i being infected by the sth infected neighbor is given as

A1=(1T1T2)s1T1,
(6)

Eq. (6) indicates that previous s − 1 infected neighbors have not changed the state of node i (not become I or SF state) before they become R state. Therefore, the probability of i being infected is

p1(R|m)=s=1mA1=T1T1+T2[1(1T1T2)m].
(7)

For case (b), node i should first become SF state, and the probability that the susceptible node i is altered by the lth infected neighbors and becomes SF can be written as

A2=(1T1T2)l1T2,
(8)

which also indicates that previous l − 1 infected neighbors have not changed the state of node i before they become R state. For the remainder ml+1 infected neighbors (including the infected neighbor who just made i go to SF state), they can infect i with probability

A3=1(1T3)ml+1.
(9)

As a result, the probability of node i first becoming SF state and then going to I state is

p2(R|m)=l=1mA2*A3=T2T1+T2[1(1T1T2)m](1T3)T2T1+T2T3[(1T3)m(1T1T2)m].
(10)

The probability p(R|m) is

p(R|m)=p1(R|m)+p2(R|m)=1(1T1T2)m(1T3)T2T1+T2T3[(1T3)m(1T1T2)m]=1T1T3+T2T3T1+T2T3(1T1T2)m(1T3)T2T1+T2T3(1T3)m.
(11)

Combining Eqs. (5) and (11), the probability of a node with degree k being infected is

p(R|k)=mp(R|m)p(m|k)=1(1θT1θT2)k(1T3)T2T1+T2T3[(1θT3)k(1θT1θT2)k]=1T1T3+T2T3T1+T2T3(1θT1θT2)k(1T3)T2T1+T2T3(1θT3)k.
(12)

Then, the EIN probability θ is the solution to the self-consistent condition

θ=kQ(k)p(R|k)=1T1T3+T2T3T1+T2T3G1(1θT1θT2)(1T3)T2T1+T2T3G1(1θT3)=f(θ).
(13)

In Eq. (13), Q(k)=(k+1)P(k+1)k is the excess degree distribution, where P(k) is the degree distribution and k is the average degree. The generating function for Q(k) is given as

G1(x)=k=0Q(k)xk.
(14)

There is a trivial solution θ = 0 in self-consistent equation (13). To have a non-trivial solution, the following condition must be met:

df(θ)dθ|θ=0=(T1+T2T3)G1(1)1,
(15)

which implies the epidemic can outbreak when

βμ+γβμ+γβ+γβFμ+β+βFkk2k.
(16)

For the prevalence of epidemic (defined as R()), we can numerically solve θ from self-consistent equation (13), then the formula of R() is

R()=kP(k)p(R|k)=1T1T3+T2T3T1+T2T3G0(1θT1θT2)(1T3)T2T1+T2T3G0(1θT3).
(17)

In Eq. (17), G0(.) is the generating function of degree distribution P(k), which is described as

G0(x)=k=0P(k)xk.
(18)

In this section, we perform an extensive set of Monte Carlo simulations to validate the theoretical predictions in Section III. Here, we carry out simulations on an Erdős-Rényi network (labeled ER network)26 with network size N = 10 000 and average degree k=10, and a configuration network generated by an uncorrelated configuration model (UCM).27 The configuration network also has N = 10 000 nodes and the degree distribution meets P(k)k3.0, whose minimal and maximal degrees are kmin = 3 and kmax=N, respectively.

Differing from the SIS (Susceptible-Infected-Susceptible) model, it is not an easy thing to determine the epidemic threshold for the SIR model owing to the non-zero value of R. In doing so, in Ref. 28, Shu et al. suggested that the variability measure

Δ=ρ2ρ2ρ
(19)

can well predict the epidemic threshold for the SIR model, where ρ denotes the prevalence of epidemic in one simulation realization.29,30 Δ can be explained as the standard deviation of the epidemic prevalence, and is a standard measure to determine critical point in equilibrium phase on magnetic system.31 In our simulations, we have taken at least 1000 independent realizations to predict the epidemic threshold. For convenience, in this study, we set recovery rate μ=1.0.

In Fig. 1, for different response rate βF, the value of R() (upper panels) and the measure Δ (lower panels) as the functions of the transmission rate β are investigated. As shown in Fig. 1, one can observe that the variability measure can well predict the epidemic threshold for our SSFIR model. As a result, in the following figures, we use this method to determine the epidemic threshold, i.e., the point where the value of Δ is the maximal. Fig. 1 also describes that, no matter γ=0.1 [see Figs. 1(a) and 1(b)] or γ=0.3 [see Figs. 1(c) and 1(d)], on the one hand, the epidemic threshold is enhanced as the response rate βF is increased. On the other hand, for a fixed value of β, Figs. 1(a) and 1(c) suggest that the prevalence of epidemic is remarkably reduced when βF is increased. The results suggest that, by introducing an additional protective state, SF, to the classical SIR model, the conclusions are quite different from the previous results which have not incorporated the impacts of human behavioral responses. The result again emphasizes the fact that the spontaneous behavioral responses of individuals to the emergent diseases have vital impacts on the epidemic dynamics. If the behavioral responses are ignored in mathematical modelling, the obtained results cannot really reflect the spreading mechanism of epidemics among human population.

FIG. 1.

On ER networks, the prevalence of epidemic R() (upper panels) and the variability measure Δ (lower panels) as functions of the transmission rate β for different values of βF and γ. Note that, here βF is a rate rather probability, as a result, whose value may be larger than one. (a)-(b) γ=0.1; (c)-(d) γ=0.3. The pink lines are given to demonstrate that the peak value of Δ corresponds to the epidemic threshold βc.

FIG. 1.

On ER networks, the prevalence of epidemic R() (upper panels) and the variability measure Δ (lower panels) as functions of the transmission rate β for different values of βF and γ. Note that, here βF is a rate rather probability, as a result, whose value may be larger than one. (a)-(b) γ=0.1; (c)-(d) γ=0.3. The pink lines are given to demonstrate that the peak value of Δ corresponds to the epidemic threshold βc.

Close modal

We then compare the theoretical results with the Monte Carlo simulations on ER network in Figs. 2 and 3. Since the degree distribution of an Erdős-Rényi network is P(k)=ekkkk!, the generating functions meet

G0(x)=G1(x)=ek(x1).
(20)

According to Ineq. (15), the epidemic threshold βc for ER network is determined by the following equation:

βcμ+γβcμ+γβc+γβFμ+βc+βF=1k.
(21)

Moreover, substituting Eq. (20) into Eqs. (13) and (17), the prevalence of epidemic R() can be easily solved.

FIG. 2.

Comparison between the Monte Carlo based simulations and the theoretical predictions for R() on ER networks. The simulation results are denoted by symbols and the theoretical predictions are denoted by the corresponding lines. The theoretical results are obtained by substituting Eq. (20) into Eqs. (13) and (17), and the value of θ is numerically solved from Eq. (13), then R() is got from Eq. (17). (a) γ=0.1; (b) γ=0.3.

FIG. 2.

Comparison between the Monte Carlo based simulations and the theoretical predictions for R() on ER networks. The simulation results are denoted by symbols and the theoretical predictions are denoted by the corresponding lines. The theoretical results are obtained by substituting Eq. (20) into Eqs. (13) and (17), and the value of θ is numerically solved from Eq. (13), then R() is got from Eq. (17). (a) γ=0.1; (b) γ=0.3.

Close modal
FIG. 3.

Comparison between the Monte Carlo based simulations and the theoretical predictions for the epidemic threshold βc on ER networks. The theoretical predictions denoted by lines are obtained from Eq. (21), and the simulation results denoted by symbols are the points where the values Δ are maximal. (a) βc as a function of the response rate βF for different values of γ; (b) βc as a function of the discount factor γ for different values of βF.

FIG. 3.

Comparison between the Monte Carlo based simulations and the theoretical predictions for the epidemic threshold βc on ER networks. The theoretical predictions denoted by lines are obtained from Eq. (21), and the simulation results denoted by symbols are the points where the values Δ are maximal. (a) βc as a function of the response rate βF for different values of γ; (b) βc as a function of the discount factor γ for different values of βF.

Close modal

The comparison of R() between the simulations and the theoretical result is plotted in Fig. 2, which indicates that the numerical simulation and the theoretical result are in good agreement. Meanwhile, the epidemic threshold for βc obtained from Eq. (21) and from numerical method (i.e., the point where Δ is maximal) is compared in Fig. 3, which also indicates that the epidemic threshold predicated by our method is remarkable agreement with numerical simulations. The result in Fig. 3 also suggests that the epidemic threshold βc is increased as the value of γ is decreased. Importantly, the reduction is more efficient when the response rate βF is larger.

Real complex networked systems often possess certain degree of skewness in their degree distributions, typically represented by some scale-free topology. We thus check our model on UCM network with degree distribution P(k)k3.0 to illustrate that our theory can generalize to the networks with heterogenous degree distribution and in the absence of degree-to-degree correlation.

As shown in Figs. 4 and 5, one can see that the analytical results are in good agreement with the numerics. They also indicate that, since increasing βF can induce more susceptible individuals go to SF state and reducing γ can lower the risk of SF nodes, as a result, both of them can lower the prevalence of epidemic R() and increase the epidemic threshold βc.

FIG. 4.

Comparison between the Monte Carlo based simulations and the theoretical predictions for R() on UCM networks. The simulation results are denoted by symbols and the theoretical predictions are denoted by the corresponding lines. The theoretical results are obtained by substituting a fixed degree distribution P(k) into Eqs. (17) and (18), and then R() can be solved from Eqs. (17) and (18) as described in Fig. 2. (a) γ=0.1; (b) γ=0.3.

FIG. 4.

Comparison between the Monte Carlo based simulations and the theoretical predictions for R() on UCM networks. The simulation results are denoted by symbols and the theoretical predictions are denoted by the corresponding lines. The theoretical results are obtained by substituting a fixed degree distribution P(k) into Eqs. (17) and (18), and then R() can be solved from Eqs. (17) and (18) as described in Fig. 2. (a) γ=0.1; (b) γ=0.3.

Close modal
FIG. 5.

Comparison between the Monte Carlo based simulations and the theoretical predictions for the epidemic threshold βc on UCM networks. The theoretical predictions denoted by lines are obtained from Ineq. (16), and the simulation results denoted by symbols are the points where the values Δ are maximal. (a) βc as a function of the response rate βF for different values of γ; (b) βc as a function of the discount factor γ for different values of βF.

FIG. 5.

Comparison between the Monte Carlo based simulations and the theoretical predictions for the epidemic threshold βc on UCM networks. The theoretical predictions denoted by lines are obtained from Ineq. (16), and the simulation results denoted by symbols are the points where the values Δ are maximal. (a) βc as a function of the response rate βF for different values of γ; (b) βc as a function of the discount factor γ for different values of βF.

Close modal

One possible way to describe our proposed SSFIR model is the mean field method, which can be written as

dSk(t)dt=βkSkΘβFkSkΘ,
(22)
dSkF(t)dt=βFkSkΘβγkSkFΘ,
(23)
dIk(t)dt=kΘβ(Sk+γSkF)Ik,
(24)
dRk(t)dt=Ik,
(25)

where the factor Θ(t)=kP(k|k)Ik(t) represents the probability that any given link points to an infected node. In the absence of any degree correlations Θ(t)=1kkkP(k)Ik(t).32 

Based on the above differential equations, we can obtain that the epidemic threshold βc=kk2 (detailed derivation is given in Section VII), which means that the epidemic threshold for our model is not related to the response rate βF or the discount factor γ, and which is the same to the epidemic threshold of classical SIR model. The simulation results based on the Eqs. (22)–(25) in Fig. 6 also indicate that, based on mean field method, the epidemic threshold is not altered by different values of βF or γ. However, our previous simulation results based on Monte Carlo method have suggested that the epidemic threshold βc is increased when βF is increased or γ is reduced. That is to say, the conclusion obtained by mean-field method is wrong.

FIG. 6.

Based on the differential equations (22)–(25), the prevalence of epidemic R() as a function of β is presented. (a) γ=0.1; (b) γ=0.3.

FIG. 6.

Based on the differential equations (22)–(25), the prevalence of epidemic R() as a function of β is presented. (a) γ=0.1; (b) γ=0.3.

Close modal

Now let us explain why the mean field method gives a wrong result. It is known that, near the epidemic threshold, the fraction of infected nodes (label ρI) is very small. When using the mean field method, the dynamic correlation is neglected, the probability of S node becoming SF is proportional to O(ρI) since the average fraction of infected nodes among neighborhood equals to O(ρI). Similarly, the probability of SF node becoming I node is also proportional to O(ρI). As a result, the probability of SSFI is proportional to O(ρI2), which leads to the effect of the SF is ignored and the epidemic threshold obtained by the mean field method is not related to the value of βF or γ. In fact, when an S node becomes an SF node there must exist at least one infected nodes among the neighborhood of the S node. More importantly, these infected neighbors may exist for a certain time interval, so the probability of SSFI is not proportional to O(ρI2). However, as deduced in Eq. (9), the dynamics correlation near the epidemic threshold is considered in our above analysis, which can accurately predict the epidemic threshold.

To summarize, we have proposed an SSFIR epidemiological model in complex networks, in which the probability of susceptible individuals becoming SF state is proportional to the number of infected neighbors, to reflect the fact that individuals are more likely to take protective measures when they find their neighbors are infected. By using theoretical analysis as well as numerical simulations, we found that the prevalence of epidemic is effectively reduced and the epidemic threshold is remarkably increased when the response rate βF is increased or the discount factor γ is reduced. Moreover, we have demonstrated that the mean field based analysis provides a wrong result: the epidemic threshold is not related to the response rate βF or discount factor γ. The reason is that, near the epidemic threshold, the probability of SSFI is a second order infinitesimal since the mean field method ignores the dynamic correlation, which makes the effect of SF state to be ignored.

With the development of technology, information induced awareness or self-protective behaviors can not only diffuse through the contact networks where the diseases spread but also can fast diffuse through many different channels, such as the word of mouth, news media, online social networks, and so on. In view of this, recent well-studied multiplex network theory may be an ideal framework to mimic the interplay of information or related awareness and the epidemic dynamics.33–35 Thus, how to generalize our model to multiplex networks and provide theoretical analysis to the model is a challenge in our further work.

This work was funded by the National Natural Science Foundation of China (Grant Nos. 61473001, 91324002, and 11331009).

When assuming Sk(0)1, then from Eq. (22), we have32 

Sk(t)=ek(β+βF)ϕ(t),
(A1)

where ϕ(t)=0tΘ(t)dt=1kkkP(k)Rk(t).

Substituting Eq. (A1) into Eq. (23), one has

dSkFdt=βFkΘ(t)ek(β+βF)ϕ(t)βγkSkFΘ(t).
(A2)

By using the variation of constants method, SkF(t) is solved as

SkF(t)=βFβ+βFβγeβγkϕ(t)(1ek(β+βFβγ)ϕ(t)).
(A3)

Then

Θ(t)=dϕ(t)dt=1ϕ(t)1kΣkP(k)ek(β+βF)ϕ(t)1kΣkP(k)βFβ+βFβγeβγkϕ(t)(1ek(β+βFβγ)ϕ(t)).
(A4)

By letting ϕ=limtϕ(t), and with the fact that limtdϕ(t)dt=0 and I(t)=0 when the epidemics end, a self-consistent equation can be got from Eq. (A4)

ϕ=11kΣkP(k)ek(β+βF)ϕ1kΣkP(k)βFβ+βFβγeβγkϕ(1ek(β+βFβγ)ϕ)=g(ϕ).
(A5)

The value ϕ=0 is always a solution. In order to have a non-zero solution, the condition

df(ϕ)dϕ|ϕ=0=(β+βF)k2kβFk2k=βk2k1
(A6)

should be satisfied, which means that the epidemic threshold βc=kk2.

1.
M. E. J.
Newman
,
SIAM Rev.
45
,
167
(
2003
).
2.
M. E. J.
Newman
,
Networks: An introduction
(
Oxford University Press
,
2010
).
3.
M.
Barthélemy
,
A.
Barrat
,
R.
Pastor-Satorras
, and
A.
Vespignani
,
J. Theor. Biol.
235
,
275
(
2005
).
4.
C.
Castellano
and
R.
Pastor-Satorras
,
Phys. Rev. Lett.
105
,
218701
(
2010
).
5.
M. E. J.
Newman
,
Phys. Rev. E
66
,
016128
(
2002
).
6.
P.
Holme
and
T.
Takaguchi
,
Phys. Rev. E
91
,
042811
(
2015
).
7.
R.
Pastor-Satorras
and
A.
Vespignani
,
Phys. Rev. Lett.
86
,
3200
(
2001
).
9.
R.
Pastor-Satorras
and
A.
Vespignani
,
Phys. Rev. E
65
,
036104
(
2002
).
10.
R.
Cohen
,
S.
Havlin
, and
D.
Ben-Avraham
,
Phys. Rev. Lett.
91
,
247901
(
2003
).
11.
H.-F.
Zhang
,
Z.
Yang
,
Z.-X.
Wu
,
B.-H.
Wang
, and
T.
Zhou
,
Sci. Rep.
3
,
3292
(
2013
).
12.
Z.
Ruan
,
M.
Tang
, and
Z.
Liu
,
Phys. Rev. E
86
,
036117
(
2012
).
13.
C. T.
Bauch
,
A. P.
Galvani
, and
D. J.
Earn
,
Proc. Natl. Acad. Sci. U. S. A.
100
,
10564
(
2003
).
14.
C. T.
Bauch
and
D. J.
Earn
,
Proc. Natl. Acad. Sci. U. S. A.
101
,
13391
(
2004
).
15.
L.
Wang
,
Y.
Zhang
,
T.
Huang
, and
X.
Li
,
Phys. Rev. E
86
,
032901
(
2012
).
16.
S.
Funk
,
M.
Salathé
, and
V. A.
Jansen
,
J. R. Soc. Interface
7
,
1247
(
2010
).
17.
X.-T.
Liu
,
Z.-X.
Wu
, and
L.
Zhang
,
Phys. Rev. E
86
,
051132
(
2012
).
18.
S.
Funk
,
E.
Gilad
,
C.
Watkins
, and
V. A.
Jansen
,
Proc. Natl. Acad. Sci.
106
,
6872
(
2009
).
19.
F. D.
Sahneh
,
F. N.
Chowdhury
, and
C. M.
Scoglio
,
Sci. Rep.
2
,
632
(
2012
).
20.
S.
Meloni
,
N.
Perra
,
A.
Arenas
,
S.
Gómez
,
Y.
Moreno
, and
A.
Vespignani
,
Sci. Rep.
1
,
62
(
2011
).
21.
Q.
Wu
,
X.
Fu
,
M.
Small
, and
X.-J.
Xu
,
Chaos
22
,
013101
(
2012
).
22.
H.-F.
Zhang
,
J.-R.
Xie
,
M.
Tang
, and
Y.-C.
Lai
,
Chaos
24
,
043106
(
2014
).
23.
N.
Perra
,
D.
Balcan
,
B.
Gonçalves
, and
A.
Vespignani
,
PLoS One
6
,
e23084
(
2011
).
24.
L.
Hébert-Dufresne
,
O.
Patterson-Lomba
,
G. M.
Goerg
, and
B. M.
Althouse
,
Phys. Rev. Lett.
110
,
108103
(
2013
).
25.
M. E. J.
Newman
and
C. R.
Ferrario
,
PLoS ONE
8
,
e71321
(
2013
).
26.
P.
Erdős
and
A.
Rényi
,
Publ. Math. Inst. Hung. Acad. Sci.
5
,
17
(
1960
).
27.
M. E. J.
Newman
,
S. H.
Strogatz
, and
D. J.
Watts
,
Phys. Rev. E
64
,
026118
(
2001
).
28.
P.
Shu
,
W.
Wang
, and
M.
Tang
, “
Simulated identification of epidemic threshold on finite-size networks
,”
Chaos
25
,
063104
(
2015
).
29.
P.
Crepey
,
F. P.
Alvarez
, and
M.
Barthélemy
,
Phys. Rev. E
73
,
046131
(
2006
).
30.
P.
Shu
,
M.
Tang
,
K.
Gong
, and
Y.
Liu
,
Chaos
22
,
043124
(
2012
).
31.
S. C.
Ferreira
,
R. S.
Ferreira
,
C.
Castellano
, and
R.
Pastor-Satorras
,
Phys. Rev. E
84
,
066102
(
2011
).
32.
Y.
Moreno
,
R.
Pastor-Satorras
, and
A.
Vespignani
,
European Phys. J. B
26
,
521
(
2002
).
33.
C.
Granell
,
S.
Gomez
, and
A.
Arenas
,
Phys. Rev. Lett.
111
,
128701
(
2013
).
34.
W.
Wang
,
M.
Tang
,
H.
Yang
,
Y.
Do
,
Y.-C.
Lai
, and
G.
Lee
,
Sci. Rep.
4
,
5097
(
2014
).
35.
S.
Boccaletti
,
G.
Bianconi
,
R.
Criado
,
C.
Del Genio
,
J.
Gómez-Gardeñes
,
M.
Romance
,
I.
Sendiña-Nadal
,
Z.
Wang
, and
M.
Zanin
,
Phys. Rep.
544
,
1
(
2014
).