Open Access+Editor's Choice • The following article is Open access

Optimal vaccination program for two infectious diseases with cross-immunity spreading on a network

Yang Ye, Qingpeng Zhang, Zhidong Cao and Daniel Dajun Zeng

Published 31 March 2021 • Copyright © 2021 The author(s)
Europhysics Letters, Volume 133, Number 4Citation Yang Ye et al 2021 EPL 133 46001DOI 10.1209/0295-5075/133/46001

Download Article PDF

Article metrics

1210 Total downloads
0 Video abstract views

Dates

Received 23 December 2020
Accepted 25 January 2021
Published 31 March 2021

Journal RSS

Abstract

There are often multiple diseases with cross-immunity competing for vaccination resources. Here we investigate the optimal vaccination program in a two-layer Susceptible-Infected-Removed (SIR) model, where two diseases with cross-immunity spread in the same population, and vaccines for both diseases are available. We identify three scenarios of the optimal vaccination program, which prevents the outbreaks of both diseases at the minimum cost. We analytically derive a criterion to specify the optimal program based on the costs for different vaccines.

Export citation and abstractBibTeX RIS

Previous article in issue

Next article in issue

Published by the EPLA under the terms of the Creative Commons Attribution 4.0 International License (CC-BY). Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

Introduction

Vaccination is an effective method of preventing epidemics like COVID-19, diphtheria, influenza, and measles [1–3]. Mathematical models have been commonly used to examine the effect of vaccination on epidemic control and prevention [4–6]. Vaccines help the human body's natural defense systems to develop antibodies to pathogens. Generally, antibodies to one pathogen do not provide protection against another pathogen. However, there is increasing evidence of cross- immunity between diseases, where exposure to one pathogen may help protect against other pathogens [7–11]. Cross-immunity is a common outcome when some diseases caused by related pathogens spread in the same population. There is rich previous epidemiological research on the interactions of pathogens with cross-immunity from the perspectives of age, gender, and population structure [12–19]. Some physical literature provides analytic calculations of the expected final epidemic sizes of multiple diseases coexisting in the same population [20,21]. It is needed to take into account the effect of cross-immunity while developing disease prevention strategies. Some studies gave simulation and optimization frameworks to identify the effective resource allocation to control several competitive diseases in the susceptible-infected-susceptible (SIS) model [22,23]. However, the identification of the optimal vaccination programs to control competitive diseases with cross-immunity in the susceptible-infected-recovered (SIR) model is under-researched and critically needed, particularly when the world is expecting to face the challenge of allocating COVID-19 vaccines in the coming years.

In this paper, we assume that there are two vaccines for two SIR-type diseases with varying cost. A vaccination program specifies the fractions of individuals that should be vaccinated for each vaccine. We define that the optimal vaccination program can control the outbreak of both diseases at the minimum cost. The microscopic Markov chain approach (MMCA) is adopted to derive the conditions guaranteeing that both diseases can be contained in the population. An optimization framework is proposed to analytically derive the optimal vaccination program. We identify three scenarios of the optimal vaccination program, and derive a criterion to specify the optimal program based on the costs for different vaccines.

Results

Here, we consider a two-layer network with the same topological structure, which represents the transmission channels for both diseases, $i=1,2$ . Two diseases propagate through different layers. Nodes represent individuals, which are susceptible (S_i), infected (I_i), or recovered (R_i) for disease i. A susceptible node S_i can be infected by an infected node I_i with an infection rate $\beta_i$ . An infected node I_i has a transition rate $\gamma_i$ to become recovered and immune to disease i. The susceptibility to disease $j (j\neq i)$ of individuals recovered from disease i is reduced by a factor $\alpha_j\in(0,1)$ (i.e., the actual infection rate for disease j becomes $\alpha_j\beta_j$ ). $1-\alpha_j$ quantifies the level of cross-immunity. In the extreme case, $\alpha_j=0$ corresponds to complete cross-immunity, and $\alpha_j=1$ corresponds to a simpler problem without cross-immunity. An individual can be simultaneously infected by both diseases. Before the occurrence of the first node infected with disease i, a fraction of $v_i\in[0,1]$ nodes are provided with vaccines for disease i. We assume that vaccine-induced immunity is equivalent to the natural immunity obtained from actual infection. Thus, vaccines for disease i directly transfer the node state from S_i to R_i without causing illness. Nodes vaccinated with vaccines for disease i are immune to disease i and also have a susceptibility to disease $j (j\neq i)$ reduced by a factor $\alpha_j$ . After vaccination, a small fraction of unvaccinated nodes become initially infected for both diseases.

According to this scheme, every node can be in nine different states at each time: S₁ S₂, S₁ I₂, S₁ R₂, I₁ S₂, I₁ I₂, I₁ R₂, R₁ S₂, R₁ I₂, and R₁ R₂. Note that the character order does not affect the state (e.g., S₁ S₂ is equivalent to S₂ S₁). Every node n has a certain probability of being in one of the nine states at time t denoted by $p^{S_1S_2}_n(t)$ , $p^{S_1I_2}_n(t)$ , $p^{S_1R_2}_n(t)$ , $p^{I_1S_2}_n(t)$ , $p^{I_1I_2}_n(t)$ , $p^{I_1R_2}_n(t)$ , $p^{R_1S_2}_n(t)$ , $p^{R_1I_2}_n(t)$ , and $p^{R_1R_2}_n(t)$ . $p^{S_1S_2}_n(t)+p^{S_1I_2}_n(t)+p^{S_1R_2}_n(t)+p^{I_1S_2}_n(t)+p^{I_1I_2}_n(t)+p^{I_1R_2}_n(t)+p^{R_1S_2}_n(t)+p^{R_1I_2}_n(t)+p^{R_1R_2}_n(t)=1$ . For disease i at time t, node n has a certain probability of being in one of three states S_i, I_i, or R_i, denoted by

where $i\neq j$ . $p_n^{R_i}(0)=v_i$ due to the effect of vaccination. If node n is recovered from disease j, the probability of not being infected for disease $i (i\neq j)$ is

where a_mn is the adjacency matrix on both layers. If node n is not recovered from disease j, the probability is

We develop the microscopic Markov chain approach to analyze the state transitions for each node [24,25]. Denote $X_i^n(t)$ as the health status of node n for disease i at time t, $X_i^n(t)\in\{S_i, I_i, R_i\}$ , then the probability that the state for node n transits from $X_1^n(t)X_2^n(t)$ to $X_1^n(t+1)X_2^n(t+1)$ can be represented as follows:

where $\mathcal{P}[X_i^n(t)\rightarrow X_i^n(t+1)]$ represents the probability that node n's health status for disease i transits from $X_i^n(t)$ to $X_i^n(t+1)$ . We present the transitions between health status for a single disease in fig. 1. Therefore, the state transitions for node n can be represented as

which is equivalent to

Fig. 1: Refer to the following caption and surrounding text. — **Fig. 1:** Illustration of the transitions between health status of node n for disease i if (a) and (b) $X_j^n(t)\neq R_j$ at time t. The definitions of $q_{n,i}^P(t)$ and $q_{n,i}^U(t)$ are given in eq. (2) and eq. (3), respectively.
Download figure:
Standard image

$X_j^n(t)=R_j$ — **Fig. 1:** Illustration of the transitions between health status of node n for disease i if (a) and (b) $X_j^n(t)\neq R_j$ at time t. The definitions of $q_{n,i}^P(t)$ and $q_{n,i}^U(t)$ are given in eq. (2) and eq. (3), respectively.
Download figure:
Standard image

The fractions of susceptible, infected, and recovered nodes at time t for disease i are denoted by $\rho^{S_i}(t)$ , $\rho^{I_i}(t)$ , and $\rho^{R_i}(t)$ , respectively,

Here N is the number of nodes on the network. The fraction of uninfected nodes at the stationary state (when t → ∞) is $\rho^{R_i}(0)+\rho^{S_i}(\infty)=v_i+\rho^{S_i}(\infty)$ for disease i. We perform extensive Monte Carlo (MC) simulations to validate the MMCA results obtained by eqs. (6). We adopt the synchronous updating method in MC simulations ¹ . All nodes update their states simultaneously in each time step, which is set as 1. We consider a network of 2000 nodes with a Poisson degree distribution, where the mean degree is 20. The initial fraction of infected nodes is $\frac{1}{N}=0.0005$ for diseases i if $v_i\neq 1$ . Each point in fig. 2 is obtained by averaging 1000 MC simulations. The average R² is 0.96, indicating a high agreement between MMCA and MC simulation results.

Fig. 2: Refer to the following caption and surrounding text. — **Fig. 2:** Comparison of the fraction of uninfected nodes at the stationary state [ $v_i+\rho^{S_i}(\infty)$ , ] using MC simulations and MMCA. We consider a network of 2000 nodes with a Poisson degree distribution, where the mean degree is 20. Results are obtained by averaging 1000 MC simulations. (a) while changing the values of v₁ and (b) while changing the values of v₂. Parameter values $\beta_1=0.4$ , $\gamma_1=0.8$ , $\alpha_1=0.1$ , $\beta_2=0.5$ , $\gamma_2=0.6$ , and $\alpha_2=0.1$ .
Download figure:
Standard image

Next, we explore the conditions preventing the outbreaks of both diseases. Specifically, the outbreak of disease i can be prevented when

According to eq. (1) and eq. (6),

Thus,

If $0\leq v_i<1$ , $p_n^{I_i}(0)=\frac{1}{N}$ , $p_n^{S_i}(0)=1-v_i-\frac{1}{N}$ , and $p_n^{R_i}(0)=0$ . Thus, $p_n^{S_iS_j}=(1-v_i-\frac{1}{N})(1-v_j-\frac{1}{N})$ , $p_n^{S_iI_j}=(1-v_i-\frac{1}{N})\frac{1}{N}$ , and $p_n^{S_iR_j}=(1-v_i-\frac{1}{N})v_j$ . If N is large enough, $p_n^{I_i}(0)=\frac{1}{N}\rightarrow0$ , consequently, from eq. (2) and eq. (3), we obtain

where k_n is the degree of node n. Thus,

where $\bar{k}$ is the mean degree of the network. To ensure the failure of outbreaks for both diseases, the right-hand side of eq. (13) should be less or equal to 0, because it is slightly larger than the left-hand side. Therefore, the following conditions should be met:

where $R_{0,1}=\frac{\beta_1\bar{k}}{\gamma_1}$ and $R_{0,2}=\frac{\beta_2\bar{k}}{\gamma_2}$ . R_{0, 1} and R_{0, 2} are the basic reproduction numbers for disease 1 and disease 2, respectively. If $v_i=1$ , $\rho_n^{I_i}(0)=\rho_n^{S_i}(0)=\rho_n^{R_i}(0)=0$ , disease i can never spread out and eqs. (14) still hold. Here R_{v, i} is the expected number of infections caused by an individual infected with disease i when the fractions of vaccinated individuals are v₁ and v₂ at the beginning. Note that this corresponds to $R_{v,i}=R_{0,i}$ for $v_1=v_2=0$ . According to eqs. (14), in order to prevent outbreaks of both diseases, the following conditions should be met:

From eqs. (15), we observe that fewer vaccines for disease i are needed to prevent the outbreak of disease i if more nodes are vaccinated for disease $j (j\neq i)$ . When $v_j=1$ (all nodes are vaccinated for disease j), $v_i \geq 1-\frac{1}{R_{0,i}\alpha_i}$ , indicating that the minimum fraction of nodes that should be vaccinated for disease i to prevent the outbreak of disease i is $1-\frac{1}{R_{0,i}\alpha_i}$ . We define this value as

In the following, we consider the common scenario when R_{0, i}>1 and $\tilde{v}_i\geq0$ . The results can be easily extended to other scenarios.

Now we plot the full $v_2-v_1$ phase diagram of the outbreak conditions for both diseases in fig. 3. We adopt the same two-layer network setting as fig. 2. As illustrated in fig. 3, only points of (v₁, v₂) in region II can prevent the outbreaks of both diseases. The intersection $(\hat{v}_1, \hat{v}_2)$ satisfies that

Denote $\Delta =[\frac{1-\alpha_2}{R_{0,1}}-\frac{1-\alpha_1}{R_{0,2}}]^2+\alpha_1^2\alpha_2^2+2\alpha_1\alpha_2[\frac{1-\alpha_2}{R_{0,1}}+\frac{1-\alpha_1}{R_{0,2}}]$ ,

Fig. 3: Refer to the following caption and surrounding text. — **Fig. 3:** Full phase diagram for the same multiplex in fig. 2. Region I: both diseases can outbreak; Region II: both diseases cannot outbreak; Region III: disease 1 cannot outbreak, while disease 2 can outbreak; Region IV: disease 2 cannot outbreak, while disease 1 can outbreak. The transition lines are computed from eqs. (15). The intersection $(\hat{v}_1, \hat{v}_2)$ is computed from eqs. (18). Parameters are as in fig. 2.
Download figure:
Standard image

$v_2-v_1$ — **Fig. 3:** Full phase diagram for the same multiplex in fig. 2. Region I: both diseases can outbreak; Region II: both diseases cannot outbreak; Region III: disease 1 cannot outbreak, while disease 2 can outbreak; Region IV: disease 2 cannot outbreak, while disease 1 can outbreak. The transition lines are computed from eqs. (15). The intersection $(\hat{v}_1, \hat{v}_2)$ is computed from eqs. (18). Parameters are as in fig. 2.
Download figure:
Standard image

Next, we aim to identify the minimum vaccination cost that can prevent the outbreaks of both diseases. Denote C₁ and C₂ as the costs of vaccination per capita for disease 1 and disease 2, respectively. C₁ and C₂ are positive. The total vaccination cost is $V=C_1v_1N+C_2v_2N$ . Then, we formulate the following optimization problem:

Since the slopes for v₁ and v₂ are both positive, the optimal vaccination program $(v_1^*$ , $v_2^*)$ is on the boundary of region II in fig. 3.

On the solid line (the boundary between region II and region IV),

and $v_2\in[\hat{v}_2, 1]$ , thus,

Therefore, V is a concave function of v₂. The optimal vaccination program $(v_1^*$ , $v_2^*)$ is $(\hat{v}_1, \hat{v}_2)$ or $(\tilde{v}_1, 1)$ . Similarly, on the dotted line (the boundary between region II and region III), the optimal vaccination program , is $(\hat{v}_1, \hat{v}_2)$ or $(1, \tilde{v}_2)$ . Overall, there are three potential scenarios of the optimal vaccination program: (1, $\tilde{v}_2$ ), ( $\tilde{v}_1$ , 1), and $(\hat{v}_1, \hat{v}_2)$ .

Comparing the values of V on these points, we find that the value of $(v_1^*, v_2^*)$ is decided by the following four nonnegative parameters:

If $L_1L_2\geq L_3L_4$ , then

thus,

If L₁ L₂<L₃ L₄, then

thus,

According to the criterion above, we can analytically derive the optimal vaccination program by substituting $(\hat{v}_1, \hat{v}_2)$ and $\tilde{v}_i$ to eq. (24) and eq. (25). In fig. 4, the optimal vaccination programs derived by analytical results are compared with the grid search with respect to different vaccine cost ratios $\frac{C_1}{C_2}$ for different diseases. The agreement is high. In fig. 4(a), $R_{0,1}=10$ , $R_{0,2}=8$ , $\alpha_1=0.4$ , $\alpha_2=0.2$ , so $L_1L_2=0.0607>L_3L_4=0.0221$ . As a result, there are three scenarios according to the vaccine cost ratio. In fig. 4(b), $R_{0,1}=10$ , $R_{0,2}=5$ , $\alpha_1=0.1$ , $\alpha_2=0.2$ , so $L_1L_2=0.0943<L_3L_4=0.3201$ . As a result, there are two scenarios according to the vaccine cost ratio.

Fig. 4: Refer to the following caption and surrounding text. — **Fig. 4:** Comparison of the optimal vaccination programs derived by analytical results and grid search with respect to different vaccine cost ratios $\frac{C_1}{C_2}$ for different diseases. Here, analytical results are calculated based on eq. (24) and eq. (25). Parameters values (a) $R_{0,1}=10$ , $R_{0,2}=8$ , $\alpha_1=0.4$ , $\alpha_2=0.2$ ; (b) $R_{0,1}=10$ , $R_{0,2}=5$ , $\alpha_1=0.1$ , $\alpha_2=0.2$ .
Download figure:
Standard image

**Fig. 4:** Comparison of the optimal vaccination programs derived by analytical results and grid search with respect to different vaccine cost ratios $\frac{C_1}{C_2}$ for different diseases. Here, analytical results are calculated based on eq. (24) and eq. (25). Parameters values (a) $R_{0,1}=10$ , $R_{0,2}=8$ , $\alpha_1=0.4$ , $\alpha_2=0.2$ ; (b) $R_{0,1}=10$ , $R_{0,2}=5$ , $\alpha_1=0.1$ , $\alpha_2=0.2$ .
Download figure:
Standard image

Conclusion

In summary, we investigate the optimal vaccination program that can prevent the outbreaks of two SIR-type diseases with cross-immunity at the minimum cost. We adopt MMCA to develop an optimization framework to identify the optimal vaccination programs in various epidemiological parameter settings. We analytically derive the optimal solutions and validate the results with MC simulations and grid search. The proposed model can be easily extended to the scenario where multiple diseases spread on the same population by quantifying the cross immunity policy between diseases. The proposed optimization framework serves as a general scheme for policy making, and can consider more than two diseases by introducing new cost functions and constraints for additional diseases. Our model can also be adapted in a sequential decision-making scheme by allowing the budget to be split into multiple sequential expenses. In such a scenario, we can incorporate the dynamic programming [26] algorithms to derive the optimal sequential decisions that maximize the final output (i.e., containing multiple diseases with minimal cost). Overall, this study provides clues to design effective and efficient vaccination programs where the cross-immunity between multiple diseases exists. In particular, the world is facing critical challenges in confronting the coexistence of both the novel coronavirus (COVID-19) pandemic and other major infectious diseases [27]. With limited resources, it is important to inform the model-driven vaccination programs that can contain multiple epidemics efficiently.

The model has limitations. First, we assume that vaccines for a specific disease provide full immunity against the disease. In practice, the efficacy of a vaccine may vary. In addition, vaccine-induced immunity may not be permanent. Besides, we do not consider the evolution of viruses in the model. In the real-world application, model calibrations are needed if the aforementioned factors are significant. Second, due to the length limit, our model does not consider the individual heterogeneity in contact patterns, susceptibility to diseases, and responses to vaccines. Incorporating the individual heterogeneity could be interesting future work. Third, we adopt a linear function to quantify the cost of the vaccine program. Quantifying the cost of vaccination program involves sophisticated economic theories and empirical data, which are not readily available yet. The cost function shall be calibrated while applying our model to public health practice.

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China (Grants No. 71972164, No. 72025404, and No. 71621002) and the Collaborative Research Fund of the Research Grants Council of Hong Kong (Grants No. C1143-20G).

Footnotes

1
The source codes are available on Github: https://github.com/jianan0099/VAC_.

Please wait… references are loading.