Local immunization program for susceptible-infected-recovered network epidemic model

The immunization can be regarded as a response to the seriousness of epidemic spreading through voluntary vaccination¹ or interventions.² Novel insights into immunization programs can be gained from the study of complex networks^3,4 from two main perspectives⁵ among others:

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

• 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 immunization^13–15 or other modelling approaches. For example, Shaban et al.¹⁶ 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 Wang¹⁷ proposed a strategy to immunize the neighbors of an infected node. Ruan et al.¹⁰ studied an information-driven vaccination program and found that strengthening the information diffusion can reduce the final vaccination fraction. Wang et al.¹⁸ investigated the interplay between information spreading and disease dynamics in an information-driven vaccination program. Jo and Baek¹⁹ and Fu et al.²⁰ evaluated efficiencies of immunization schemes in the susceptible-infected-recovered-susceptible (SIRS) and susceptible-infected-susceptible (SIS) networks, respectively. Zhang et al.²¹ investigated the impact of subsidy policies on vaccination decisions under the voluntary vaccination. More recently, Takaguchi et al.²² 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).²³ It becomes of interest to evaluate the cost-effectiveness of static as well as dynamic immunization schemes for this realistic situation. A recent paper²³ 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 nodes¹⁷ 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:¹⁹ 

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 = k₀, $k0+1, ⋯,$ k_c, where k₀ and k_c are the minimal and maximal degrees. Assuming I_k(0) ≃ 0 and R_k(0) = 0 for each k, then the probability Θ(t) of a randomly selected node connecting to an infected individual can be formulated as^27,28

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

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

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

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,⁸ 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

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 I_k(0) ≃ 0 and R_k(0) = 0 for each k, the occurrence of an epidemic outbreak depends on the stability of the disease free equilibrium of the network model.³² Hence, we consider the system near the zero solution (I_k(t) = 0 for each k). Then, $E[w2(s)]≃β(1−δ)kΘ$ and the evolution of the infected class³³ can be given by the linearized model of (2)

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

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).³² 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.³⁴ It is also worthy to remark that the same disease outbreak threshold (4) can also be obtained by the branching process theory.³⁵ 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.³⁶ 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.¹⁴ 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;³⁸ 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.

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 resource³ 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 G₁ formed by the nodes in Ω can be defined, and the remaining nodes together with their edges form a subnetwork G₂. Since there exist links between nodes in G₁ and G₂, the whole network G can be regarded as an interdependent network³⁹ where the mean-field approach is still feasible. Taking the difference between interdependent networks and a single network, we call the subnetwork G_i, 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 approach⁴⁰).

The random immunization means that a fraction f of all nodes is randomly selected to be immunized,³ from which, the set Ω is determined. Based on the nodes in set Ω, blocks G₁ and G₂ 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 G_i (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

Here, the respective infection probabilities in G₁ and G₂ are

and

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 E₀ gives

with $Θ=⟨k⟩−1∑k(k−1)P(k)[Ik(1)+Ik(2)]$⁠.

At the early stage of an epidemic, we have

while

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

with

where $I2M×2M$ is an identity matrix and A is a M × M positive matrix with entries $Akk′=k(k′−1)P(k′)/⟨k⟩$ and M = k_c − k₀ + 1.

Using the property of block matrix, one can obtain

Using this equality, it is easy to obtain the maximal eigenvalue of J₁(E₀) as

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

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.³ To evaluate the efficacy of this program in this study, we choose nodes with large degrees to be vaccinated, that is, the node set $Ω={v∈W: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:

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

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 E₀

The corresponding Jacobian matrix J₂(E₀) at E₀ becomes

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

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

with

Therefore, the epidemic threshold for the targeted immunization case is

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

with $f̃=∑k≥Kk(k−1)P(k)∑kk(k−1)P(k)$⁠. In the targeted immunization program, we have $f=∑k≥KP(k)$⁠, which represents total vaccine coverage in the whole network. Next, we are going to show that $f̃>f$⁠, which is equivalent to

Since k(k − 1) is an increasing function of k, there exists a threshold value m such that $k(k−1)≤∑j=k0kcj(j−1)P(j)$ when k ≤ m and $k(k−1)≥∑j=k0kcj(j−1)P(j)$ when k ≥ m, indicating that the function F increases first and then decreases across the threshold value. On the other hand, F(k₀) = 0 and F(k_c) > 0, and therefore, we get F(K) > 0 for all K > k₀, which proves $f̃$ > f. Therefore, $τct>τcr$ for the same f value with $f=∑k≥KP(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 model⁴¹ with degree exponent α = 2.7. The network structure is set with size N = 2000, the minimal degree k₀ = 3, and the maximal degree $kc=N≃44.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 I_max for different parameters, which has been shown to be an effective index to measure the epidemic threshold of the model with infinite absorbing states.³³ An alternative approach is based on the variability measure suggested by Shu et al.⁴² 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.

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

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

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.²³ 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.

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,⁴³ the branching process approximation,¹⁶ 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).⁴⁵ 

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-Satorras⁴⁶ 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 networks³⁹ or the community networks.⁴⁷ 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:

and

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

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.²³