Responsive immunization and intervention for infectious diseases in social networks

The spread of an epidemic disease (e.g., tuberculosis (TB),¹ severe acute respiratory syndrome (SARS),² Asian-influenza,³ swine-origin influenza A (H1N1)^4,5) in a population can be studied by using dynamical system approaches. The susceptible-infected-susceptible (SIS) model and the susceptible-infected-recovered/removed (SIR) model are two widely considered models. The theory of complex network can provide an analytic framework for heterogenous contact patterns of a population.⁶ The heterogeneous contact reflects the property that the node degree k, or, the number of contacts with other individuals for a given individual, is not uniform. The frequently used network models are Erdös-Rényi (ER) random graphes,⁷ Strogatz-Watts (SW) small-world networks,⁸ and Barabási-Albert (BA) scale-free networks.⁹ There is a sharp contrast in their degree distributions P(k). The study of epidemic spreading on complex networks shows the existence of a high correlation between the condition of epidemic outbreak and the degree distribution.^10–14 

How to find optimal immunization strategies to minimize the risk of epidemic outbreaks on complex networks have been widely studied.^14–21 A number of basic immunization strategies have been proposed and investigated, such as the random immunization,¹⁴ the targeted immunization¹⁴ and the acquaintance immunization.¹⁵ These proposed immunization strategies, however, are built on a major premise that vaccination or immunization is compulsory and have not considered the willingness or desire of individuals. Given some social factors, such as, religious belief and human rights, thus, the immunization behavior is not a compulsory behavior but decided by individuals themselves. In this situation, whether to vaccinate or not is related to the risk of being infected by the infectious disease.^18,20 As a kind of individual awareness,^22,23 risk assessment of infection should be closely related to an individual's local information (the state of an individual's neighborhood). Hence, individual vaccination behavior depends on its local information. For the convenience, we call such dynamic immunization as the responsive immunization,²¹ which is also referred as the information-driven vaccination¹⁹ or the information dependent vaccination.²⁴ 

In our work²² and further work,²⁵ the local information affects individual susceptibility and can change the epidemic threshold. An interesting problem is: can the responsive immunization based on the local information affect the epidemic threshold? Intuitively, the answer seems to be “yes.” However, we will find that it is not the case for the responsive immunization. In this paper, we mainly focus on the influence of the responsive immunization on the “epidemic threshold.” Herein, the “epidemic threshold” means the critical value that can discriminate the dynamical behaviors of the system.²⁶ 

The rest of this paper is organized as follows: in Sec. II, an SIS epidemic model with responsive immunization strategy is proposed and the theoretical analysis on the thresholds of the model is given based on the microscopic Markov-chain method. Numerical simulations are also presented to verify the theoretical results. In Sec. III, the external intervention in the responsive immunization is further investigated. At last, conclusions and some discussions are summarized in Sec. IV.

We use the SIS model to investigate the effect of responsive immunization. The SIS model is chosen for the following two reasons: (1) it is widely applicable and may be adapted for some epidemic diseases such as meningitis and gonorrhea;²⁷ (2) it is also relatively simple and approximated to the early stage of the epidemic outbreak. We also assume an epidemic spreads along the static network G with size N, which is completely determined by its adjacency matrix A where the entries $aij=aji=1$ if there is a link between nodes i and j, otherwise $aij=aji=0$⁠. All nodes of G are enumerated with index $i=1,2,⋯,N$⁠.

In this model with the responsive immunization considered in the SIS model, each individual may stay in one of three states: S-susceptible, I-infected, and V-vaccinated. During a time step, a susceptible individual may get infected at an average rate β per unit time if it is contacted by one infected individual, and meanwhile may also be vaccinated and then removed (due to the responsive immunization) at rate $πvac(i)$⁠. An infected individual may recover and become susceptible again at rate γ per unit time. It is assumed that all these events are independent.

Intuitively, the immunization rate $πvac$ increases with both the local information of an epidemic and the response rate of epidemic risk, δ (⁠$0≤δ≤1$⁠). Suppose node i with degree k_i has $kinf$ (⁠$kinf≥0$⁠) infected neighbors, then the immunization rate of node i is given by

When $δ=0$⁠, the model is just the standard SIS epidemic model.¹⁰ Similar to the previous work, we define the effective spreading rate $λ=β/γ$⁠. For the sake of the following analysis, we first present a lemma:

Lemma 1: Let $V*$ be a set composed of nodes in a network G and $pj,t$ be the probability of node $j∈V*$ to be infected at time t, respectively. Then the number of infected nodes in $V*$ is a stochastic variable $ξ∈[0,|V*|]$ (⁠$|V*|$ denoting the number of the elements in set $V*$⁠) and its expected value satisfies

Proof: Note that

According to the inductive method, we can complete the proof.

Previous literatures have indicated that the microscopic Markov-chain approximation (MMA) approach^28,29 is an effective method in studying the epidemic spreading in quenched networks (i.e., the adjacency matrix is unchanged in time), including unweighted networks,^26,30 weighted networks,^29,31 and even multiplex networks.³² 

Along this way, we denote the probability of node i to be infected and to be immunized at time t, $pi,t$ and $qi,t$⁠, respectively. During the time interval $[t,t+1)$⁠, the change of $pi,t$ depends on two events: the recovery from state I to state S and the infection from state S to state I. The change of $qi,t$ only depends on the vaccination from state S to state V.

Note that for each susceptible node i, there exists two exclusive events per unit time: (i) getting infected at rate $πinf(i)$⁠; (ii) becoming immunized at rate $πvac(i)$⁠. Therefore, with these notations the discrete-time epidemic network model is described by

Now, we establish the specific forms of $πinf(i)$ and $πvac(i)$⁠. Let $ζi,t$ denote the probability of node i being uninfected at time t. Considering node i may get infected through connections with each of its infected neighbors, the expression of uninfected probability reads^28,29

Here, $Ni$ denotes the neighborhood of node i. Accordingly, we have

We continue to give the expression for $πvac(i)$⁠. Following its definition, we consider the stochastic variable $ξ=kinf$ (⁠$ξ≤ki$⁠). By using Lemma 1, we take the expected value of ξ as an approximation to $kinf$ and obtain

So Eq. (1) can be rewritten as

When $δ=0$ and $qi,0=0$ for each i, model (5) reduces into the following simple form:

This is just the standard SIS model without reinfection terms.^28,29 According to the previous literature, the critical value of epidemic outbreak for the spreading rate, $λcSIS$⁠, obeys

where $Λmax(A)$ is the leading eigenvalue of the adjacency matrix $A$⁠. When the spreading rate is larger than the critical value, the epidemic disease will become endemic and persist in a population.

For the approximation model (5), we first analyze the steady state of its dynamical behavior. To this end, we would like to determine the values of $pi,t$ and $qi,t$ at the steady state, p_i and q_i. On substituting $pi,t≡pi$ and $qi,t≡qi$ into Eqs. (5), we have

As long as $δ>0$⁠, we have $(1−pi−qi)∑j∈Nipj=0$⁠. This indicates two possible cases: (i) $∑j∈Nipj=0$⁠; (ii) $1−pi−qi=0$⁠. When $∑j∈Nipj=0$⁠, we have $pj=0$ for each $j∈Ni$ and further $1−∏j∈Ni(1−βpj)=0$⁠. In either case, $(1−pi−qi)[1−∏j∈Ni(1−βpj)]≡0$⁠. This means that for each node i, $pi=0$ but the value of q_i is unknown and related to the initial conditions of the system. Therefore, the fraction of infected node in a population (denoted by I(t)) always decays to zero regardless of the spreading rate λ and the response rate δ (⁠$δ>0$⁠).

Next, we want to estimate the critical value of epidemic spreading for our model, that is, the epidemic threshold $λc$⁠. Based on the analysis of the steady state, the epidemic threshold means: if $λ≤λc$⁠, I(t) decreases to zero (no epidemic), otherwise, first increases to a maximum and then decreases to zero (an epidemic) due to the responsive immunization. Following,^33,34 the occurrence or not of an epidemic depends on the stability of the disease free equilibrium of the disease model described by (5).

Additionally, we notice that $I(t)=(1/N)∑i=1Npi,t$⁠. Therefore, in order to obtain the mathematical expression of the epidemic threshold we only consider the subsystem (5a) near the disease free equilibrium (⁠$pi=0$ and $0≤qi≤1$ for each i). At this time, a linear form of Eqs. (5a) can be written as

This uses the approximation $(1−a)(1−b)≃1−a−b$ when $a≪1,b≪1$⁠.

This system is not a closed form since we do not know q_i value corresponding to $pi=0$ for each i. However, we can approximately analyze the stability of the disease free equilibrium under the assumption (H): $pi,0≃0$ and $qi,0=0$ for each i. In fact, from Eqs. (5), when $pi,0=0$ we have $pi,t=0$ for $t>0$ and then $qi,t=0$⁠. In other words, under the assumption (H) $pi,0→0$ can imply that $qi,t→0$ and further $qi→0$⁠. Hence, we consider the system near the zero solution (⁠$pi=qi=0$ for each i). At this time, model (8) becomes

Now we study the stability of the above system. Let us introduce a vector function $pt=(p1,t,p2,t,...,pN,t)T∈RN$ (the state vector of the network). Then, Eq. (9) can be given by a collective form

Here, $I$ denotes a N-dimensional identity matrix. From above, the local stability of the zero solution of system (10) can be established by

This indicates that the epidemic threshold for the spreading rate, $λc$⁠, is given by

The threshold condition is the same as that in the standard SIS model (6) and indicates that the response rate δ has no impact on the epidemic threshold. This differs significantly from our previous work without vaccination,²² where we have shown that the local information can affect the epidemic threshold. We argue that the epidemic threshold is unchanged for different δ values because: (1) the response rate does not directly affect the dynamic of infection but decrease the number of susceptible nodes; (2) at the beginning of an epidemic spreading, the vaccination fraction generated by the responsive immunization is not large enough to halt the epidemic outbreak.

It is worth stressed that the analysis of the epidemic threshold is on the basis of the assumption (H). If the assumption is not satisfied, the epidemic threshold should be related with initial conditions (i.e., $pi,0$⁠, $qi,0$⁠). In addition, considering the effectiveness of the MMA approach, we believe that the responsive immunization has no significantly affect on the epidemic threshold in the stochastic model, together with the following simulation analysis.

To test above argument, we perform Monte Carlo simulations over both BA scale-free networks⁹ with the degree distribution $P(k)∼k−3$ and ER networks with connecting probability $p=0.006$⁠.⁷ The Monte Carlo simulations are implemented in a parallel way, that is, each node's state can be updated with a certain rate in a time step. More specifically, for a susceptible node i at each time step, we generate a random number $r∈[0,1)$⁠, if $r∈[0,πinf(i))$ then node i is changed from state S to state I; else if $r∈[πinf(i),πinf(i)+πvac(i))$ then node i is changed from state S to state V; else node i is still susceptible. Here, $πinf(i)$ and $πvac(i)$ are the transition rates as stated above. One should note that, since we mainly focus on the epidemic threshold of the model, near this critical point, the epidemic just begins to prevail, so the infected neighbors are few, in this situation, we can image that the value of $πinf(i)+πvac(i)$ is smaller than 1. Although the network considered here is small (N = 2000), we also made simulations for a larger network (e.g., N = 5000) and obtained a similar observation.

We first verify the accuracy of Eq. (11). Simulations begin with a single seed initial condition. To minimize random fluctuation caused by the initial conditions, we make average over 200 realizations of different initial infected nodes.

In Fig. 1, we compare the Monte Carlo simulations and theoretical model (5) on the epidemic threshold, $λc$⁠. The Markov chain prediction shows that the epidemic threshold is indeed independent of the parameter δ, which complies with the analysis from Eq. (11). In order to further verify this formula, we compute the value of $Λmax(A)$ for the BA scale-free network used here according to the power method and obtain that $Λmax(A)≃26.325$⁠. Following Eq. (11), we have that $λc=1−Λmax(A)≃0.03799$⁠, which is a good approximation to the simulation results obtained from Eqs. (5) (Fig. 1(b)). For an ER random network, $Λmax(A)≃13.202$⁠, so $λc=1−Λmax(A)≃0.07575$⁠. This also complies with the simulation results obtained from Eq. (5) (Fig. 1(d)).

The epidemic threshold for the Monte Carlo simulation on a BA scale-free network (Fig. 1(a)) is larger than that for the Markov chain prediction based on Eqs. (5) (Fig. 1(b)). This is mainly due to the first order approximation of the mathematic model.³⁵ 

Hinted by the real situations in which the responsive immunization often takes place after the beginning of the outbreak of the epidemic (for an emerging infectious disease at least), we consider another kind of initial condition: a portion of infected seeds, e.g., 1% of the nodes are infected. Unless otherwise specified, we set the recovery rate $γ=0.5$ in the later simulations.

Herein, we would like to examine the exactness of the theoretical model (5) compared to the Monte Carlo simulation. In order to do this, we consider the maximal infection density $Imax$ and the corresponding peak time $tmax$ over different BA scale-free networks and ER random networks with the almost same mean degree, respectively. Fig. 2 shows the peak time and the maximal infection density as functions of the infection rate for different δ values, which are consistent with theoretical models as indicated by the solid lines. This figure tells us that the MMA approach is effective to model the epidemic spreading with responsive immunization, even for the prediction of the peak time.

Furthermore, Fig. 2 suggests that, though the epidemic threshold is not related to the value of δ, increasing the value of δ can lower the peak time $tmax$⁠, indicating the upward tendency of epidemic is fast controlled. Meanwhile, the maximal infection density $Imax$ also decreases with the value of δ, which means that the responsive immunization can effectively control or hinder the prevalence of epidemic.

In the final part of this section, we investigate the initial outbreak of epidemic, which means that I(t) increases at t = 0. For the discrete-time system (5), it is needed to consider the quantity

Considering the average over initial conditions, we have $pi,0=ε$ and $qi,0=0$ for each i. Then,

On plugging the above equality into Eq. (12), we obtain

From this equality, the critical value of initial epidemic outbreak obeys

Interestingly, we find that $λc0$ is inversely proportional to the mean degree of the network. When $λ>λc0$⁠, $ΔI>0$ and I(t) first increases; when $λ<λc0$⁠, we have $ΔI<0$ and the epidemic prevalence first decreases. One should note that, this condition is used to judge whether the epidemic will outbreak at initial time step, though $λ<λc0$⁠, the outbreak of epidemic is also possible to happen since $λc0>λc$ (as illustrated in the latter simulation).

In Fig. 3, we illustrate the change of $ΔI$ on a BA scale-free network with the mean degree $⟨k⟩≅12$⁠. Note that $ε=0.01$ and $γ=0.5$⁠, so we have $ΔI=−0.005+0.12β.$ This tells us that $ΔI$ is a linear function of the infection rate β, which can be seen in this figure. Also, we can derive the condition of initial epidemic outbreak: $λ>λc0⇒β>0.0417$⁠. In the inset of Fig. 3, we can see that $ΔI>0$ only if $β≥0.05$⁠, which is in accordance with the Markov chain prediction.

According to the content of the responsive immunization, each individual potentially vaccinates with a certain rate. Hence, the vaccine coverage may range from nodes with small degree to ones with large degree. In order to investigate the range of the vaccine coverage in our model, we consider the degree distribution of vaccinated nodes (called as the vaccination degree distribution)

which can be compared to the degree distribution of all nodes in the network P(k) by using numerical simulations.

In this section, we use BA scale-free networks to simulate epidemic dynamics and initially 1% of the nodes are infected. As shown in Fig. 4, the relation between two kinds of degree distribution is $F(k)∼P(k)∼k−3$⁠. The good relation should be induced by the degree uncorrelated property of BA networks and it would be interesting to consider the impact of network structure on the vaccination degree distribution and its relation with the degree distribution in other work. Nevertheless, this tells us that the vaccination possibility of the node with small degree is nearly equal to that of the node with large degree, which forces a large amount of vaccines to be required for the responsive immunization (see Fig. 5(a): c = 0). However, the situation that too many people choose vaccination potentially leads to the waste of resources. Especially when the vaccine is rare, it may cause social panic, even violence. So we should take some necessary intervention to the voluntary vaccination.

It is well known that the targeted immunization scheme in scale-free networks is very effective in controlling the epidemic outbreak.¹⁴ Inspired by this, we introduce the intervention measure based on the targeted immunization and only permit those nodes with large degree to take vaccine—when node i goes to the epidemic prevention station to vaccinate, we only allow those nodes with degree $ki≥c$ to get vaccinated. In other words, we modify the above model with responsive immunization as follows:

Here, c reflects the level of vaccination intervention. The larger the value of c, the higher the intervention level. When c = 0, there is no any vaccination intervention. As above analysis, an infectious disease dies out eventually. However, in this case the vaccination fraction is generally large. Hence, we consider the impact of the external intervention (i.e., $c>0$⁠) on the vaccination fraction.

In Fig. 5(a), we report the change of the vaccination fraction at the steady state $V=V∞(c)$ for different c values when the maximal response rate (i.e., $δ=1$⁠) is assumed. In this figure, we choose the parameter window $β∈[0,0.06]$ where for all cases become zero/vanish as $t→∞$ (Fig. 5(b)). This can allow us to clearly observe the intervention effect. As a result, V decreases significantly with increasing c value. We take the case $β=0.05$ as an illustration. In this case, $V∞(0)≅70%$⁠, while $V∞(10)≅30%$ and $V∞(30)≅1%$⁠.

For the model with no intervention, the infection density becomes zero eventually. But for the model with degree-based intervention, the change of infection density I(t) may not be like that. Actually, as shown in Fig. 5(b), when the infection rate is large enough there exists a positive value of infection density at the steady state. This determines another epidemic threshold $λc′$⁠. When $λ>λc′$⁠, $I=I∞(c)>0$⁠; when $λ<λc′$⁠, I = 0. Therefore, the intervention can lead to another larger epidemic threshold, above which the infectious disease will undergo a new outbreak and persist in the population. This is not desirable from the perspective of epidemic control and can be solved by increasing the vaccine coverage. When the vaccine is rare, we should implement other strategies, e.g., quarantine strategies.

Besides the SIS-like epidemic threshold $λc′$⁠, the system has two other critical values: (1) the threshold of initial epidemic outbreak $λc0$⁠; (2) the SIR-like epidemic threshold $λc$ above which the epidemic disease asymptomatically decays. We can approximately evaluate these quantities. Actually, after substituting the modified form of $πvac(i)$ (14) into model (1) and using the similar analysis, we obtain that $λc0=1/⟨k⟩,$ and $λc≅1/Λmax(A),$ where the usage of the approximation symbol accounts for the impact of a portion of initial infection condition. These equalities indicate that the two parameters δ and c almost have no impact on both $λc0$ and $λc$⁠.

Since $Λmax(A)$ is never smaller than the mean degree of the network, $⟨k⟩$⁠, $λc≤λc0$⁠.³⁰ What is more, $λc$ must be less than $λc′$⁠. When $λc′≠λc0$⁠, the dynamical behavior of the system under consideration is completely characterized by three critical values of spreading rate (⁠$λc,λc0$⁠, and $λc′$⁠). As an illustration, we consider the intervention model for the case c = 30 and $γ=0.5$ on a BA scale-free network. According to our simulations and computations, three critical values of infection rate $βc≃0.02,βc0≃0.04$⁠, and $βc′≃0.06$⁠. Therefore,

When the infection rate lies between $βc$ and $βc0$⁠, the epidemic prevalence first decreases, then increases, and finally drops to zero. Fig. 6 shows that the change of the infection curve I(t) is just like that. We also notice the deviation at the peak time between the Monte Carlo simulation and the Markov chain prediction. This may be related to the number of infected seeds distributed in the network initially. When this number is very small, a possible outbreak may be eradicated due to the simulation randomness. Hence, we can see that the deviation becomes smaller with increasing I(0) values as illustrated in Fig. 6(b).

As a brief summary, a modified SIS epidemic model with the responsive immunization is proposed and analyzed. Our main results are: (1) we derive the epidemic threshold of the model and find that $λc=1/Λmax(A)$⁠; (2) we establish the critical condition of initial epidemic outbreak; (3) we obtain the impact of the intervention on the epidemic dynamics and final vaccination size.

Recently, Sahneh and et al.²⁶ found two kinds of epidemic thresholds that can discriminate the asymptotical decay and the exponential decay of an epidemic disease in a susceptible-alert-infected-susceptible (SAIS) model. In our intervention model, we also find the existence of multiple critical values of spreading rate. Interestingly, the meanings of $λc$ and $λc′$ are just as stated in Ref. 26. From the viewpoint of epidemic control, it is not enough to raise the larger threshold $λc′$ because the epidemic still breaks out when the spreading rate λ satisfies $λc<λ<λc′$⁠.

In our model, the responsive immunization has no significant influence on the epidemic threshold and a large vaccination fraction is required to halt the epidemic spreading. It is necessary to take some intervention measures in the voluntary vaccination. We argue that there are two types of intervention measures according to their effects: (1) one is that the epidemic threshold becomes larger; (2) the other is that the epidemic threshold is unchanged but final vaccination size decreases. In the present paper, the degree-based intervention based on the targeted immunization belongs to the second type since it can only reduce the vaccine size. One may consider the first type of vaccination intervention, e.g., the responsive immunization on the adaptive contact network¹⁷ where the epidemic threshold may be changed due to the change of network structure.

The authors wish to thank the referees whose comments improved the quality of the paper. This research was supported by the National Natural Science Foundation of China (Nos. 61203153 and 11005001).