The immunization strategies through contact tracing on the susceptible-infected-recovered framework in social networks are modelled to evaluate the cost-effectiveness of information-based vaccination programs with particular focus on the scenario where individuals belonging to a specific set can get vaccinated due to the vaccine shortages and other economic or humanity constraints. By using the block heterogeneous mean-field approach, a series of discrete-time dynamical models is formulated and the condition for epidemic outbreaks can be established which is shown to be not only dependent on the network structure but also closely related to the immunization control parameters. Results show that increasing the immunization strength can effectively raise the epidemic threshold, which is different from the predictions obtained through the susceptible-infected-susceptible network framework, where epidemic threshold is independent of the vaccination strength. Furthermore, a significant decrease of vaccine use to control the infectious disease is observed for the local vaccination strategy, which shows the promising applications of the local immunization programs to disease control while calls for accurate local information during the process of disease outbreak.

Researchers have performed many studies about the dynamic immunization schemes on networks, which provide new insights into optimal strategies to effectively control the diseases. However, previous work mainly focused on the immunization over the whole population under the assumption that the vaccines supply is sufficient to cover the most of individuals at risk. In some special scenarios, especially for emerging diseases, the vaccine supply is small compared with the large amount of population size at risk. Therefore, it is pivotal to investigate a local immunization program that is confined to a specified subpopulation. In the current manuscript, the immunization strategies through contact tracing on the susceptible-infected-recovered framework in social networks are proposed. Using dynamic analysis and numerical simulations, the epidemic thresholds are found which are shown to be closely associated with the adjustable parameters for the vaccination programs. This result is different from the case for susceptible-infected-susceptible model framework where epidemic threshold is independent of the vaccination strength. Furthermore, a significant decrease of vaccine use to control the infectious disease, through comparison of static and dynamic immunization schemes, is observed for the local vaccination strategy. These results provide novel designs for disease control using immunization programs.

The immunization can be regarded as a response to the seriousness of epidemic spreading through voluntary vaccination1 or interventions.2 Novel insights into immunization programs can be gained from the study of complex networks3,4 from two main perspectives5 among others:

  • The static immunization which is implemented before the epidemic spreading.6 Generally speaking, there are two basic schemes following this idea, the random immunization and the targeted immunization,3 and other multiple variant strategies, such as the acquaintance immunization7 and inverse targeting immunization,6 with applications to time-varying networks8 and multiplex networks.9 

  • The dynamic immunization where the program is implemented during the epidemic outbreak. Since an individual's vaccination decision is made mainly based on the epidemic seriousness while sometimes, the novel effective and safe vaccine can only be developed and mass-produced after the emergence of infectious diseases triggered by new pathogens, the dynamic immunization is always implemented in realistic situations in most cases.1,10

Similar to individual's behavioral responses to infectious diseases,11,12 the dynamic immunization can be well-adjusted based on transmission information during the epidemic spreading. As soon as an infectious disease begins to spread in the population, an effective immunization program should be initiated which may be adjusted according to the disease prevalence and the program terminates when the disease dies out. Therefore, the epidemic information-based immunization allows us to take the advantage of the interplay between the immunization response and the epidemic spreading. Motivated by this idea, researchers have evaluated the efficacy of various dynamic immunization schemes by using network models with direct immunization13–15 or other modelling approaches. For example, Shaban et al.16 formulated real-time susceptible-infected-recovered (SIR) vaccination models for contact tracing in a network with a specific degree distribution by a branching process approximation. Nian and Wang17 proposed a strategy to immunize the neighbors of an infected node. Ruan et al.10 studied an information-driven vaccination program and found that strengthening the information diffusion can reduce the final vaccination fraction. Wang et al.18 investigated the interplay between information spreading and disease dynamics in an information-driven vaccination program. Jo and Baek19 and Fu et al.20 evaluated efficiencies of immunization schemes in the susceptible-infected-recovered-susceptible (SIRS) and susceptible-infected-susceptible (SIS) networks, respectively. Zhang et al.21 investigated the impact of subsidy policies on vaccination decisions under the voluntary vaccination. More recently, Takaguchi et al.22 proposed an immunization strategy based on observer placement, which is shown to be very efficient for disease control in the clustering networks.

Although providing novel insights into the cost-effectiveness of various immunization approaches, previous studies mainly focused on the immunization over the whole population (denoted by W in this manuscript). However, in reality, the vaccine supply is limited compared with the large population size at risk, especially for those diseases triggered by novel pathogens. Furthermore, vaccines may deliver to only a partial population due to other economic or humanity constraints. Therefore, it becomes much more realistic to consider a local immunization program that only covers a specified subpopulation (denoted by Ω in what follows).23 It becomes of interest to evaluate the cost-effectiveness of static as well as dynamic immunization schemes for this realistic situation. A recent paper23 investigated the SIS epidemic model with local immunization program and showed that the condition of epidemic outbreak is not related to the immunization strength. The SIS modelling framework is well-accepted for describing infections, such as rotavirus and gonorrhea, which do not confer long-lasting immunity. However, a SIR framework is more suitable for infections such as measles, mumps, and chickenpox where individuals recover and confer lifelong immunity, which is the fundamental framework we will extend in the current paper with network structures. Furthermore, we are also interested in the comparison of predictions obtained from the immunization programs for these two network model frameworks.

In the current manuscript, we are going to propose two kinds of hypothetical immunization models, for local and global immunization programs, applicable to the SIR epidemic networks. Rigorous and numerical analysis will illustrate the disease transmission conditions, and the results will further be compared with previous studies on vaccination programs on the SIS framework (see the  Appendix and Ref. 23). The rest of this paper is organized as follows: In Section II, we first formulate a susceptible-vaccinated-infected-recovered (SVIR) model with a global immunization program; then in Section III, we further investigate a local immunization program by theoretical analysis and simulations; finally, discussions are presented in Section IV.

In this section, we extend the SIR epidemic model to SVIR models with the consideration of a global immunization scheme. Herein, the V state represents the immunized population through vaccination. Intuitively, one possible efficient immunization strategy is to directly immunize all the susceptible nodes connected to infected nodes, refereed to as the high-risk nodes17 which are likely to be infected by their infectious neighbors in the following infection wave. These high-risk nodes can be found through contact tracing theoretically.16,24,25 However, in practice, it is not easy to locate all these high-risk nodes and there exists a discount rate. Hence, we introduce an adjustable parameter δ, denoted as the tracing rate or immunization rate, to account the efficiency of tracing and immunizing these high-risk nodes. Suppose a susceptible node will get infected by one infectious neighbor with rate β. Then a susceptible node with s infected neighbors changes its state with the following probabilities:19 

P(SS)=(1δ)s(1β)s,P(SV)=1(1δ)s:=w1(s),P(SI)=(1δ)s[1(1β)s]:=w2(s).
(1)

Following Ref. 26, we assume the network is randomly generated according to the degree distribution P(k)kα with α ∈ (2, 3] as many real networks to incorporate the heterogeneity of individuals. This assumption implies that the connectivity of nodes is uncorrelated. We denote Sk(t),Vk(t),Ik(t),Rk(t) as the relative densities of susceptible, vaccinated, infected, and recovered nodes in the population with degree k at time step t, respectively, with k = k0, k0+1,,kc, where k0 and kc are the minimal and maximal degrees. Assuming Ik(0) ≃ 0 and Rk(0) = 0 for each k, then the probability Θ(t) of a randomly selected node connecting to an infected individual can be formulated as27,28

Θ(t)k(k1)P(k)Ik(t)kkP(k)=k(k1)P(k)Ik(t)k.

Then the probability that a node of degree k has exactly s infected neighbors is given by the binomial distribution29 

Bin(k,s)=(ks)Θs(1Θ)ks.

Taking the expectation of the stochastic variable w1(s) with respect to the above defined binomial distribution gives the probability with which a susceptible node of degree k is vaccinated

E[w1(s)]=1sBin(k,s)(1δ)s=1(1δΘ)k.

Similarly, a susceptible node of degree k gets infected with probability

E[w2(s)]=sBin(k,s)(1δ)s[1(1β)s]=(1δΘ)k[1(δ+βδβ)Θ]k.

In the present paper, we employ the widely used discrete-time approach,29–31 capable of accounting the periodicity feature in daily life or day-night changes,8 in the above process of state changes during disease transmission. If we assume that an infected node recovers and simultaneously achieves the perpetual immunization to the pathogen with rate γ, then the discrete-time epidemic process can be described in a mean-field form

Sk(t+1)=Sk(t)Sk(t){E[w1(s)]+E[w2(s)]},Vk(t+1)=Vk(t)+Sk(t)E[w1(s)],Ik(t+1)=(1γ)Ik(t)+Sk(t)E[w2(s)],Rk(t+1)=Rk(t)+γIk(t).
(2)

Clearly, variables Sk(t),Vk(t),Ik(t),Rk(t) are nonnegative and satisfy Sk(t)+Vk(t)+Ik(t)+Rk(t)=1 for each k and t, which can be shown from the system (2), or by their definitions.

Under the initial conditions Ik(0) ≃ 0 and Rk(0) = 0 for each k, the occurrence of an epidemic outbreak depends on the stability of the disease free equilibrium of the network model.32 Hence, we consider the system near the zero solution (Ik(t) = 0 for each k). Then, E[w2(s)]β(1δ)kΘ and the evolution of the infected class33 can be given by the linearized model of (2)

Ik(t+1)=(1γ)Ik(t)+β(1δ)kΘ.
(3)

By analyzing the Jacobian matrix of Eqs. (3), one can find that the system has a unique eigenvalue of maximum modulus, i.e., β(1δ)k1kk(k1)P(k)+1γ, from which the epidemic threshold can be derived

τc=11δkk2k,
(4)

with k2=kk2P(k). Here, we use the rescaled infection rate τ = β/γ. Then the epidemic threshold τc determines the epidemic outbreak: if τ < τc, the total infection density I(t)=kIk(t)P(k) decreases to zero (no epidemic), otherwise, I(t) first increases to a maximum and then decreases to zero (an epidemic).32 When δ > 0, the epidemic threshold is inversely proportional to 1 − δ value.

It is interesting to observe that the threshold index τc derived here equals to the epidemic threshold for the SIR model with local information-based behavioral responses.34 It is also worthy to remark that the same disease outbreak threshold (4) can also be obtained by the branching process theory.35 For example, the authors in Ref. 16 studied the vaccination through the contact tracing with a general contact time and obtained the similar result, while τ = β/(β + γ) when the contact time follows an exponential distribution in that paper.

The development of optimal vaccine allocation strategies to control the epidemic spreading remains a central problem in public health and network security.36 Furthermore, vaccine shortages, resulting from higher-than-expected demand, interruptions in production/supply or a lack of budgets, makes it impossible to immunize almost all the nodes in the whole network and urges to design an optimal strategy minimizing the total number of vaccines or the social cost.14 The design of an optimal strategy in the consideration of this constraint is not only related to the high-risk nodes but also the nodes with other particular characteristics during the epidemic spreading.6,37 For that purpose, we introduce a local immunization program that is confined to a node set Ω, in which nodes are predefined according to special characteristics for the epidemic control, with the special reference to vaccine shortage scenario: only susceptible nodes in Ω can be vaccinated or removed while other nodes cannot. An illustrative diagram is shown in Fig. 1, where the infected node at the center is surrounded by 8 susceptible nodes and only one of two traced susceptible nodes in Ω gets vaccinated. Clearly, ΩW with two extreme cases: when Ω=, the immunization scheme through vaccination is not implemented and the model reduces to a standard SIR model;38 while the global immunization program is implemented to cover all risky nodes for Ω = W. In what follows, we only consider the local immunization scenario, that is, the above inclusions are proper.

FIG. 1.

Illustrations of the contact tracing by an infected node. The central node is an infected node with 8 susceptible neighbors, among which two are in the set Ω while the others are outside Ω. Only one, out of two Ω-nodes, is traced (and also will get vaccinated) by the central node along the contact between them (indicated by a blue and thick line).

FIG. 1.

Illustrations of the contact tracing by an infected node. The central node is an infected node with 8 susceptible neighbors, among which two are in the set Ω while the others are outside Ω. Only one, out of two Ω-nodes, is traced (and also will get vaccinated) by the central node along the contact between them (indicated by a blue and thick line).

Close modal

In the local immunization program, one key question is to determine which group of individuals should be traced and get immunized, in other words, to define the set Ω for the network. This question can be solved when we have no knowledge about the spreading resource3 by some classical static immunization strategies, such as the random immunization and targeted immunization. In order to present a comparative analysis, in this paper, we consider two kinds of Ω determined by the random immunization and targeted immunization, respectively.

When the set Ω is fixed, a subnetwork G1 formed by the nodes in Ω can be defined, and the remaining nodes together with their edges form a subnetwork G2. Since there exist links between nodes in G1 and G2, the whole network G can be regarded as an interdependent network39 where the mean-field approach is still feasible. Taking the difference between interdependent networks and a single network, we call the subnetwork Gi, i = 1, 2, as blocks of G and the mean-field approach based on the blocks is correspondingly referred to as the block heterogeneous mean-field (HMF) approach (it differs from the block variable mean-field approach40).

The random immunization means that a fraction f of all nodes is randomly selected to be immunized,3 from which, the set Ω is determined. Based on the nodes in set Ω, blocks G1 and G2 can be defined accordingly. Denote Sk(i)(t),Vk(i)(t),Ik(i)(t), and Rk(i)(t) as the relative densities of susceptible, vaccinated, infected, and recovered nodes of Gi (i = 1 or 2) in the population G with degree k at time step t, respectively. According to the discrete-time HMF approach, the dynamical model for the random immunization program is given by

Ik(1)(t+1)=(1γ)Ik(1)(t)+Sk(1)(t)E[w2(1)(s)],Ik(2)(t+1)=(1γ)Ik(2)(t)+Sk(2)(t)E[w2(2)(s)].

Here, the respective infection probabilities in G1 and G2 are

w2(1)(s)=(1δ)s[1(1β)s],

and

w2(2)(s)=1(1β)s.

Please note that the model does not include the dynamics of other variables for Sk(i)(t),Rk(i)(t) with i = 1, 2 and Vk(1) which do not appear in the system describing the evolution of infectious nodes at the initial stage of disease spread.

A similar approximation analysis as in Sec. II near the disease-free equilibrium E0 gives

E[w2(1)(s)]β(1δ)kΘandE[w2(2)(s)]βkΘ,

with Θ=k1k(k1)P(k)[Ik(1)+Ik(2)].

At the early stage of an epidemic, we have

Sk(1)(t)fandSk(2)(t)1f

while

Ik(i)(t)0,Rk(i)(t)0foreachiandVk(1)(t)0.

Let I(t)=(Ik0(1),Ik0+1(1),,Ikc(1),Ik0(2),Ik0+1(2),,Ikc(2))T, then the local stability of E0 can be established through the following linear system for infected nodes:

I(t+1)=J1(E0)I(t)

with

J1(E0)=(1γ)I2M×2M+β(f(1δ)Af(1δ)A(1f)A(1f)A),

where I2M×2M is an identity matrix and A is a M × M positive matrix with entries Akk=k(k1)P(k)/k and M = kc − k0 + 1.

Using the property of block matrix, one can obtain

det(f(1δ)AλIM×Mf(1δ)A(1f)A(1f)AλIM×M)=det((1fδ)AλIM×M0(1f)AλIM×M)=λMdet[λIM×M(1fδ)A].

Using this equality, it is easy to obtain the maximal eigenvalue of J1(E0) as

λmax(J1)=1γ+β(1fδ)λmax(A).

Therefore, the epidemic threshold for the random immunization case, τcr, is given by

τcr=11fδkk2k.
(5)

It is obvious to see that the epidemic threshold is dependent of both f and δ. In addition, the epidemic threshold becomes much more sensitive to the immunization rate δ for larger f values.

The targeted immunization has been shown very effective in controlling epidemic outbreak on scale-free networks.3 To evaluate the efficacy of this program in this study, we choose nodes with large degrees to be vaccinated, that is, the node set Ω={vW:deg(v)K}, where deg(v) denotes the degree of node v and K is a control parameter. In this scenario, the whole network can be divided into two blocks:

G1withdegreekKandG2withdegreek<K.

The dynamical equations for the targeted immunization case can be written as

Ik(1)(t+1)=(1γ)Ik(1)(t)+Sk(1)(t)E[w2(1)(s)],kK,Ik(2)(t+1)=(1γ)Ik(2)(t)+Sk(2)(t)E[w2(2)(s)],k<K.

Here, E[w2(1)(s)] and E[w2(2)(s)] are defined as before. We can obtain the linearized equations for I(t)=[Ik0(2),Ik0+1(2),...,IK(2),IK+1(1),...,Ikc(1)]T for linear stability analysis of the disease free equilibrium E0

I(t+1)=J2(E0)I(t).

The corresponding Jacobian matrix J2(E0) at E0 becomes

J2(E0)=(1γ)IM×M+βB,
(6)

where B is a M × M positive matrix with entries Bkk=k(k1)P(k)/k for k < K and Bkk=(1δ)k(k1)P(k)/k for k ≥ K.

It is easy to verify from (6) that the dominant eigenvalue

λmax(J2)=1γ+βλmax(B),

with

λmax(B)=k<Kk(k1)P(k)k+(1δ)kKk(k1)P(k)k.

Therefore, the epidemic threshold for the targeted immunization case is

τct=kk<Kk(k1)P(k)+(1δ)kKk(k1)P(k)=kk2kδkKk(k1)P(k).
(7)

It indicates that the epidemic threshold is dependent on δ and K. The epidemic threshold increases with δ while decreases with K, and the infectious disease may be controlled when δ is large or K is small enough. To compare the cost-effectiveness of the random and targeted immunization strategies, we first write Eq. (7) into the form of Eq. (5) as

τct=11f̃δkk2k
(8)

with f̃=kKk(k1)P(k)kk(k1)P(k). In the targeted immunization program, we have f=kKP(k), which represents total vaccine coverage in the whole network. Next, we are going to show that f̃>f, which is equivalent to

F(K):=k=Kkc[ k(k1)kk(k1)P(k) ]P(k)>0.

Since k(k − 1) is an increasing function of k, there exists a threshold value m such that k(k1)j=k0kcj(j1)P(j) when k ≤ m and k(k1)j=k0kcj(j1)P(j) when k ≥ m, indicating that the function F increases first and then decreases across the threshold value. On the other hand, F(k0) = 0 and F(kc) > 0, and therefore, we get F(K) > 0 for all K > k0, which proves f̃ > f. Therefore, τct>τcr for the same f value with f=kKP(k) for the targeted immunization case, which implies that the epidemic threshold becomes greater for the target immunization program than that for the random immunization program with the same vaccine coverage used. Therefore, one can conclude that the targeted immunization is more efficient than the random immunization.

To verify the above theoretical analysis, we perform Monte Carlo simulations over scale-free networks generated from the standard configuration model41 with degree exponent α = 2.7. The network structure is set with size N = 2000, the minimal degree k0 = 3, and the maximal degree kc=N44.72. The recovery rate γ is set to be 1.0. All simulations are implemented by a parallel updating strategy in which the actual disease states of each node and its neighbors at each time step are considered. We start with a single initial infectious seed and all simulation results are obtained by taking averages of 20 random network configurations and 50 independent initial conditions for each network realization.

Fig. 2(a) illustrates the epidemic prevalence R (i.e., the fraction of recovered nodes at the end of the epidemic wave) as a function of infection rate β. This figure also shows the existence of an epidemic threshold for different immunization rates. In order to examine the validation of the theoretical results to the Monte Carlo simulation, we also consider the maximal infection density Imax for different parameters, which has been shown to be an effective index to measure the epidemic threshold of the model with infinite absorbing states.33 An alternative approach is based on the variability measure suggested by Shu et al.42 As illustrated in Fig. 2(b), the simulation results agree with the theoretical threshold conditions obtained in Eq. (5). Similar conclusion can be made for the targeted immunization in Fig. 3. Furthermore, increasing δ value can always raise the epidemic threshold τc regardless of the random or targeted immunization case. This result is significantly different from the SVIS network model based on an SIS framework (see detailed analysis about SVIS model in the  Appendix), where the parameter δ does not play a role for τc.

FIG. 2.

Effect of the immunization rate on the epidemic threshold and prevalence for the random immunization case when f = 0.5: (a) final recovered size R versus the infection rate β for different values of δ; (b) contour of Imax in the (δβ) parameter plane, where the white line indicates the theoretically predicted curve determined in Eq. (5), and the light gray region corresponds to the parameter region with zero prevalence.

FIG. 2.

Effect of the immunization rate on the epidemic threshold and prevalence for the random immunization case when f = 0.5: (a) final recovered size R versus the infection rate β for different values of δ; (b) contour of Imax in the (δβ) parameter plane, where the white line indicates the theoretically predicted curve determined in Eq. (5), and the light gray region corresponds to the parameter region with zero prevalence.

Close modal
FIG. 3.

Effect of the immunization rate on the epidemic threshold and prevalence for the targeted immunization case when K = 10: (a) final recovered size R versus the infection rate β for different values of δ; (b) contour of Imax in the (δβ) parameter plane, where the white line indicates the theoretically predicted curve determined in Eq. (7), and the light gray region to the parameter region with zero prevalence.

FIG. 3.

Effect of the immunization rate on the epidemic threshold and prevalence for the targeted immunization case when K = 10: (a) final recovered size R versus the infection rate β for different values of δ; (b) contour of Imax in the (δβ) parameter plane, where the white line indicates the theoretically predicted curve determined in Eq. (7), and the light gray region to the parameter region with zero prevalence.

Close modal

We further investigate the impact of dynamic immunization on the final vaccine size V by using the immunization efficiency Q. When there is no infection in the network (i.e., at the steady state), we can define ΩX as

ΩX={iΩ|state{i}=X},

where X is the node state, which may be S, I, R, or V. Notice that there exist infection-induced immunization nodes in Ω and therefore, ΩR ∪ ΩS ∪ ΩV = Ω. Hence, the immunization efficiency for the SIR model, a function of variables δ, β, f, and K, can be expressed as

Q(δ,β,f)=Ω#SΩ#Ω#R,

where the symbol #ΩX denotes the number of the elements in set ΩX.

Fig. 4 clearly shows that the immunization efficiency index Q is strongly correlated with the infection rate β and the immunization rate δ for different predefined sets Ω, representing the random/targeted immunization strategies used. For each fixed δ value, Q increases as β decreases. However, the monotonicity of Q, as a function of δ is much more complicated, as shown in Fig. 4(a) for the random immunization case. Generally speaking, Q is an increasing function of δ (Figs. 4(b)–4(d)). However, in Fig. 4(a), when the infection rate β is relatively small, say β = 0.2, the immunization efficiency is negatively correlated to the immunization strength, as highlighted by three blue lines. It is due to the dual effect of the increased immunization strength.23 Although increasing δ enlarges vaccination coverage, it also halts the spreading of an epidemic with a small infection rate across hub nodes and hence decreases propensity for vaccination. The relationship between τc and δ is illustrated by dashed lines in each panel of Fig. 4. Almost all of these curves lie in the red region, showing that the immunization efficiency should be very high to control the disease. It is noticed that at the case δ = 1, the expression of τc is reduced to be the same as the corresponding static immunization in Ref. 38. However, the spreading patterns between dynamic and static immunization are not the same, as many susceptible nodes are not vaccinated in Ω for the dynamic immunization.

FIG. 4.

The contour plot of the immunization efficiency, where the horizontal coordinate is the infection rate β and the vertical coordinate is the immunization strength δ. Panels (a) and (b) illustrate the random immunization case for f = 0.2 and 0.8, respectively, while panels (c) and (d) show the targeted immunization case for K = 6 and 12, respectively. The dashed lines in each panel indicate the epidemic thresholds by theoretical predictions.

FIG. 4.

The contour plot of the immunization efficiency, where the horizontal coordinate is the infection rate β and the vertical coordinate is the immunization strength δ. Panels (a) and (b) illustrate the random immunization case for f = 0.2 and 0.8, respectively, while panels (c) and (d) show the targeted immunization case for K = 6 and 12, respectively. The dashed lines in each panel indicate the epidemic thresholds by theoretical predictions.

Close modal

The study of the network theory enables us to analyze the role of each node or node set in the epidemic spreading and get novel insights into the transmission dynamics. As we know, the SIR-like epidemic network model can be analyzed by various approaches, such as the percolation theory,43 the branching process approximation,16 and the effective degree approaches.32,44 Recently, the heterogeneous mean-field approach poses a good tool to analyze complex disease dynamics due to its simple and deterministic formulation.18,28 In this manuscript, we formulated real-time immunization models with the discrete-time HMF approach, where susceptible nodes can get immunized by contact tracing from infected nodes. Considering the real situations of vaccine shortages such that the number of vaccines cannot cover the whole population, we propose a local immunization program that can only immunize a given node set Ω in the whole population, which can be defined as a geological region of a city or a social group of a population or other groups sharing some characteristics. The epidemic thresholds for different (local versus global) vaccination scenarios against infectious diseases are obtained from stability analysis, based on which the effectiveness of a vaccination program can be evaluated. The predicted thresholds are validated through numerical simulations. Our result suggests that the local immunization program can greatly improve the efficiency of static immunization, requiring a smaller amount of vaccines to effectively control disease spread. However, the efficacy of vaccination programs not only depends on immunization rate but also on the choice of individual group to immunize. Therefore, it remains pivotal to extend the approach in this manuscript to other local immunization strategies, with different targeted vaccination groups to get an optimal strategy for disease control. This may contribute toward the optimal strategy of vaccine allocation for emerging infectious diseases such as influenza A (H1N1).45 

In the local immunization program with a SIR framework, we find that the immunization rate δ can greatly affect the epidemic threshold, which distinguishes from the prediction based on the SIS spreading mechanism where δ does not play a role in the threshold. This adds one more difference between the SIS and SIR network models, as revealed by Castellano and Pastor-Satorras46 that the threshold of generic epidemic models is vanishing for an SIS model, while it is finite for the SIR model on quenched scale-rich networks (i.e., α > 3).

In the present work, we only consider the same immunization rate δ for each node in the immunization set Ω. The same approach remains valid for a general case with multiple immunization sets with different δ values. Another interesting exploration may be the study of local immunization program in the interdependent networks39 or the community networks.47 These realistic issues suggest good topics for further research.

Q.W. would like to thank National Natural Science Foundation of China for its support under No. 61203153. Y.L. was supported in part by NSFC (11301442) and RGC (PolyU 253004/14P). We thank the referees for their supportive comments which greatly improve the whole manuscript.

In this appendix, we extend the SIS epidemic model with the consideration of dynamic immunization programs. For comparative purposes, the immunization and disease transmission parameters are set to be the same as those in the SVIR model. We denote Sk(i)(t),Vk(i)(t), and Ik(i)(t) as the susceptible, vaccinated, infected densities among nodes inside (i = 1) and outside (i = 2) of the tracing group Ω with degree k at time step t, respectively. Then the following discrete-time SVIS model can be formulated for two different groups by using a similar approach as that in the main part:

Sk(1)(t+1)=Sk(1)(t){ 1E[w1(1)(s)]E[w2(1)(s)] }+γIk(1)(t),Vk(1)(t+1)=Vk(1)(t)+Sk(1)(t)E[w1(1)(s)],Ik(1)(t+1)=(1γ)Ik(1)(t)+Sk(1)(t)E[w2(1)(s)],
(A1)

and

Sk(2)(t+1)=Sk(2)(t){ 1E[w2(2)(s)] }+γIk(2)(t),Ik(2)(t+1)=(1γ)Ik(2)(t)+Sk(2)(t)E[w2(2)(s)].
(A2)

In this model, w1(1)(s)=1(1δ)s,w2(1)(s)=(1δ)s[1(1β)s], and w2(2)(s)=1(1β)s.

In order to obtain an approximate condition for epidemic outbreak, we set Sk(i)(t)=Sk(i),Ik(i)(t)=Ik(i) for i = 1, 2 and Vk(1)(t)=Vk(1) in Eqs. (A1) and (A2) where Sk(i),Ik(i), and Vk(1) are constant values at the equilibrium state and furthermore, assume Ik(1)=0 such that there is no infection in the set Ω due to immunization. Then the infectious proportion Ik(2) satisfies

γIk(2)=[Nk(2)Ik(2)][1(1βΘ)k],

where Nk(2)=Sk(2)+Ik(2) is a fixed value for each k and Θ=kkP(k)Ik(2)kkP(k). The positive solution for Ik(2) of the above algebraic equation exists if and only if τ:=βγ>kkP(k)kNk(2)k2P(k). This inequality gives the epidemic threshold for the SVIS model, which is independent of the immunization strength δ. A similar result was claimed in the previous work.23 

1.
T. C.
Reluga
and
A. P.
Galvani
,
Math. Biosci.
230
,
67
(
2011
).
2.
Task Force on Community Preventive Services
,
Am. J. Prev. Med.
18
,
92
(
2000
).
3.
R.
Pastor-Satorras
and
A.
Vespignani
,
Phys. Rev. E
65
,
036104
(
2002
).
4.
M. E. J.
Newman
,
Networks: An Introduction
(
Oxford University Press
,
Oxford
,
2010
).
5.
J.
Goldenberg
,
Y.
Shavitt
,
E.
Shir
, and
S.
Solomon
,
Nat. Phys.
1
,
184
(
2005
).
6.
M.
Schneider
,
T.
Mihalijev
, and
H. J.
Herrmann
,
Europhys. Lett.
98
,
46002
(
2012
).
7.
R.
Cohen
,
S.
Havlin
, and
D.
Ben Avraham
,
Phys. Rev. Lett.
91
,
247901
(
2003
).
8.
B. A.
Prakash
,
H.
Tong
,
N.
Valler
,
M.
Faloutsos
, and
C.
Faloutsos
, in
Proceedings of the 2010 European Conference on Machine Learning and Knowledge Discovery in Databases: Part III, ECML PKDD'10
(
2010
).
9.
D. W.
Zhao
,
L. H.
Wang
,
S. D.
Li
,
Z.
Wang
,
L.
Wang
, and
B.
Gao
,
PLoS One
9
,
e112018
(
2014
).
10.
Z.
Ruan
,
M.
Tang
, and
Z. H.
Liu
,
Phys. Rev. E
86
,
036117
(
2012
).
11.
S.
Funk
,
M.
Salathé
, and
V. A. A.
Jansen
,
J. R. Soc. Interface
7
,
1247
(
2010
).
12.
S. F.
Chen
and
Q. C.
Wu
,
J. Jiangxi Normal Univ. (Nat. Sci. Ed.)
39
,
531
(
2015
).
13.
X. L.
Peng
,
X. J.
Xu
,
X. C.
Fu
, and
T.
Zhou
,
Phys. Rev. E
87
,
022813
(
2013
).
14.
L.
Chen
and
J.
Sun
,
Physica A
410
,
196
(
2014
).
15.
D. Q.
Shi
,
L.
Ke
,
J. G.
Zhou
, and
R. Z.
Yu
,
J. Jiangxi Normal Univ. (Nat. Sci. Ed.)
37
,
637
(
2013
).
16.
N.
Shanban
,
M.
Andersson
,
A.
Svensson
, and
T.
Britton
,
Math. Biosci.
216
,
1
(
2008
).
17.
F. Z.
Nian
and
X. Y.
Wang
,
J. Theor. Biol.
264
,
77
(
2010
).
18.
W.
Wang
,
M.
Tang
,
H.
Yang
,
Y.
Do
,
Y. C.
Lai
, and
G.
Lee
,
Sci. Rep.
4
,
5097
(
2014
).
19.
H. T. M. H. H.
Jo
and
S. K.
Baek
,
Physica A
361
,
534
(
2006
).
20.
X. C.
Fu
,
M.
Small
,
D. M.
Walker
, and
H. F.
Zhang
,
Phys. Rev. E
77
,
036113
(
2008
).
21.
H. F.
Zhang
,
Z. X.
Wu
,
X. K.
Xu
,
M.
Small
,
L.
Wang
, and
B. H.
Wang
,
Phys. Rev. E
88
,
012813
(
2013
).
22.
T.
Takaguchi
,
T.
Hasegawa
, and
Y.
Yoshida
,
Phys. Rev. E
90
,
012807
(
2014
).
23.
Q. C.
Wu
,
X.
Fu
,
Z.
Jin
, and
M.
Small
,
Physica A
419
,
566
(
2015
).
24.
L. S.
Tsimring
and
R.
Huerta
,
Physica A
325
,
33
(
2003
).
25.
I. Z.
Kiss
,
D. M.
Green
, and
R. R.
kao
,
Proc. R. Soc. B
1570
,
1407
(
2005
).
26.
R.
Pastor-Satorras
and
A.
Vespignani
,
Phys. Rev. Lett.
86
,
3200
(
2001
).
27.
M.
Barthélemy
,
A.
Barrat
,
R.
Pastor-Satorras
, and
A.
Vespignani
,
J. Theor. Biol.
235
,
275
(
2005
).
28.
I. Z.
Kiss
,
D. M.
Green
, and
R. R.
Kao
,
Math. Biosci.
203
,
124
(
2006
).
29.
30.
S.
Gómez
,
A.
Arenas
,
J.
Borge-Holthoefer
,
S.
Meloni
, and
Y.
Moreno
,
Europhys. Lett.
89
,
38009
(
2010
).
31.
Y.
Shang
,
Int. J. Biomath.
6
,
1350007
(
2013
).
32.
J.
Lindquist
,
J. L.
Ma
,
P.
van den Driessche
, and
F. H.
Willeboordse
,
J. Math. Biol.
62
,
143
(
2011
).
33.
Q. C.
Wu
,
H. F.
Zhang
, and
G. H.
Zeng
,
Chaos
24
,
023108
(
2014
).
34.
H. F.
Zhang
,
J. R.
Xie
,
M.
Tang
, and
Y. C.
Lai
,
Chaos
24
,
043106
(
2014
).
35.
P.
Jagers
,
Branching Processes with Biological Applications
(
Wiley-Interscience
,
London
,
1975
).
36.
A. N.
Hill
and
I. M.
Longini
, Jr.
,
Math. Biosci.
181
,
85
(
2003
).
37.
L.
Hébert-Dufresne
,
A.
Allard
,
J. G.
Young
, and
L. J.
Dubé
,
Sci. Rep.
3
,
2171
(
2013
).
38.
N.
Madar
,
T.
Kalisky
,
R.
Cohen
,
D.
Ben Avraham
, and
S.
Havlin
,
Eur. Phys. J. B
38
,
269
(
2004
).
39.
M.
Dickison
,
S.
Havlin
, and
H. E.
Stanley
,
Phys. Rev. E
85
,
066109
(
2012
).
40.
S. Y.
Liu
,
N.
Perra
,
M.
Karsai
, and
A.
Vespignani
,
Phys. Rev. Lett.
112
,
118702
(
2014
).
41.
M. E. J.
Newman
,
S. H.
Strogatz
, and
D. J.
Watts
,
Phys. Rev. E
64
,
026118
(
2001
).
42.
P.
Shu
,
W.
Wang
, and
M.
Tang
,
Chaos
25
,
063104
(
2015
).
43.
M. E. J.
Newman
,
Phys. Rev. E
66
,
016128
(
2002
).
44.
C. R.
Cai
,
Z. X.
Wu
, and
J. Y.
Guan
,
Phys. Rev. E
90
,
052803
(
2014
).
45.
M.
Baguelin
,
A. J. V.
Hoek
,
M.
Jit
,
S.
Flasche
,
P. J.
White
, and
W. J.
Edmunds
,
Vaccine
28
,
2370
(
2010
).
46.
C.
Castellano
and
R.
Pastor-Satorras
,
Phys. Rev. Lett.
105
,
218701
(
2010
).
47.
Z. H.
Liu
and
B.
Hu
,
Europhys. Lett.
72
,
315
(
2005
).