Interplay between the local information based behavioral responses and the epidemic spreading in complex networks

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.¹⁸ 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,¹⁹ 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.²⁰ 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, S^F, 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 (SS^FIR) 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 SS^FIR model²³ 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 S^F 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,⁵ our results show that the introduction of S^F 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 S^F 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 S^F, is introduced into the original SIR model, we use SS^FIR 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 S^F state with a response rate $βF≥0$⁠. The transmission rate for the S^F nodes is lower; thus, we assume the transmission rate for S^F nodes is $γβ$⁠, where $0≤γ<1$ is a discount factor.

The SS^FIR model is described by the four following reactions and the associated rates:

Note that the SS^FIR model returns to the SIR model once $βF=0$⁠, and the S^F 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, S^FI, and SR edge are β, β_F, and μ, respectively. As a result, the probabilities of an SI edge becoming an II and S^FI edge are given as $T1=ββ+βF+μ$ and $T2=βFβ+βF+μ$⁠, respectively. Similarly, since the transition rate of an S^FI edge becoming an II and S^FR during a sufficiently small time interval $Δt$ are $γβ$ and μ, the probability of an S^FI edge becoming an II edge is $T3=γβγβ+μ$⁠.²⁴ 

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 large²⁵). 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

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 S^F 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

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

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

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 $m−l+1$ infected neighbors (including the infected neighbor who just made i go to S^F state), they can infect i with probability

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

The probability $p(R|m)$ is

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

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

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

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

which implies the epidemic can outbreak when

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

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

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)²⁶ with network size N = 10 000 and average degree $〈k〉=10$⁠, and a configuration network generated by an uncorrelated configuration model (UCM).²⁷ The configuration network also has N = 10 000 nodes and the degree distribution meets $P(k)∼k−3.0$⁠, whose minimal and maximal degrees are k_min = 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

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.³¹ 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 SS^FIR 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, S^F, 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.

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)=e−〈k〉〈k〉kk!$⁠, the generating functions meet

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

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

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)∼k−3.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 S^F state and reducing γ can lower the risk of S^F nodes, as a result, both of them can lower the prevalence of epidemic $R(∞)$ and increase the epidemic threshold β_c.

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

where the factor $Θ(t)=∑k′P(k′|k)Ik′(t)$ represents the probability that any given link points to an infected node. In the absence of any degree correlations $Θ(t)=1〈k〉∑kkP(k)Ik(t)$⁠.³² 

Based on the above differential equations, we can obtain that the epidemic threshold $βc=〈k〉〈k2〉$ (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.

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 S^F is proportional to $O(ρI)$ since the average fraction of infected nodes among neighborhood equals to $O(ρI)$⁠. Similarly, the probability of S^F node becoming I node is also proportional to $O(ρI)$⁠. As a result, the probability of $S→SF→I$ is proportional to $O(ρI2)$⁠, which leads to the effect of the S^F 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 S^F 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 $S→SF→I$ 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 SS^FIR epidemiological model in complex networks, in which the probability of susceptible individuals becoming S^F 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 $S→SF→I$ is a second order infinitesimal since the mean field method ignores the dynamic correlation, which makes the effect of S^F 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 have³² 

where $ϕ(t)=∫0tΘ(t′)dt′=1〈k〉∑kkP(k)Rk(t)$⁠.

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

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

Then

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

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

should be satisfied, which means that the epidemic threshold $βc=〈k〉〈k2〉$⁠.