Optimal vaccination program for two infectious diseases with cross immunity

Vaccination is an effective method of preventing epidemics like COVID-19, diphtheria, influenza, and measles Andre et al. (2008); Bonanni (1999); Le et al. (2020). Mathematical models have been commonly used to examine the effect of vaccination on epidemic control and prevention Wang et al. (2016); Dykman et al. (2008); Khanjanianpak et al. (2020). Vaccines help the human body’s natural defense systems to develop antibodies to pathogens. Generally, antibodies to one pathogen do not provide protection to another pathogen. However, there is increasing evidence of cross immunity between diseases, where exposure to one pathogen may help protect against other pathogens. Wu et al. (2020); Schultz-Cherry (2015); Balmer and Tanner (2011); Laurie et al. (2015, 2018). 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 pathogens with cross immunity from the perspectives of age, gender, and population structure Andreasen et al. (1997); Castillo-Chavez et al. (1989, 1996); Li et al. (2003); Webb et al. (2013); Leventhal et al. (2015); Darabi Sahneh and Scoglio (2014); Azimi-Tafreshi (2016). Some physical literature provides analytic calculations of the expected final epidemic sizes of multiple diseases coexisting in the same population Karrer and Newman (2011); Newman (2005). 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 Watkins et al. (2016); Chen and Preciado (2014). 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 letter, 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.

Here, we consider a two-layer network with the same topological structure, which represents the transmission channels for both diseases, $i=1,2$, which 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 also immune to disease $i$ and 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_{1}S_{2}$, $S_{1}I_{2}$, $S_{1}R_{2}$, $I_{1}S_{2}$, $I_{1}I_{2}$, $I_{1}R_{2}$, $R_{1}S_{2}$, $R_{1}I_{2}$, and $R_{1}R_{2}$. Note that the character order does not affect the state (e.g., $S_{1}S_{2}$ is equivalent to $S_{2}S_{1}$). Every node $n$ has a certain probability of being in one of the nine states at time $t$ denoted by $p^{S_{1}S_{2}}_{n}(t)$, $p^{S_{1}I_{2}}_{n}(t)$, $p^{S_{1}R_{2}}_{n}(t)$, $p^{I_{1}S_{2}}_{n}(t)$, $p^{I_{1}I_{2}}_{n}(t)$, $p^{I_{1}R_{2}}_{n}(t)$, $p^{R_{1}S_{2}}_{n}(t)$, $p^{R_{1}I_{2}}_{n}(t)$, and $p^{R_{1}R_{2}}_{n}(t)$. $p^{S_{1}S_{2}}_{n}(t)+p^{S_{1}I_{2}}_{n}(t)+p^{S_{1}R_{2}}_{n}(t)+p^{I_{1}S_{2}}_{n}(t)+p^{I_{1}I_{2}}_{n}(t)+p^{I_{1}R_{2}}_{n}(t)+p^{R_{1}S_{2}}_{n}(t)+p^{R_{1}I_{2}}_{n}(t)+p^{R_{1}R_{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

	$\displaystyle p_{n}^{S_{i}}(t)$	$\displaystyle=p_{n}^{S_{i}S_{j}}(t)+p_{n}^{S_{i}I_{j}}(t)+p_{n}^{S_{i}R_{j}}(t),$		(1)
	$\displaystyle p_{n}^{I_{i}}(t)$	$\displaystyle=p_{n}^{I_{i}S_{j}}(t)+p_{n}^{I_{i}I_{j}}(t)+p_{n}^{I_{i}R_{j}}(t),$
	$\displaystyle p_{n}^{R_{i}}(t)$	$\displaystyle=p_{n}^{R_{i}S_{j}}(t)+p_{n}^{R_{i}I_{j}}(t)+p_{n}^{R_{i}R_{j}}(t),$

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

$$q_{n,i}^{P}(t)=\prod\limits_{m}(1-a_{mn}p_{m}^{I_{i}}(t)\alpha_{i}\beta_{i}),$$

(2)

where $a_{mn}$ is the adjacency matrix on both layers. If node $n$ is not recovered from disease $j$, the probability is

$$q_{n,i}^{U}(t)=\prod\limits_{m}(1-a_{mn}p_{m}^{I_{i}}(t)\beta_{i}).$$

(3)

We develop the microscopic Markov chain approach to analyze the state transitions for each node $n$ as

$$\begin{split}p_{n}^{S_{1}S_{2}}(t+1)&=q_{n,1}^{U}(t)q_{n,2}^{U}(t)p_{n}^{S_{1}S_{2}}(t),\\ p_{n}^{S_{1}I_{2}}(t+1)&=q_{n,1}^{U}(t)[1-q_{n,2}^{U}(t)]p_{n}^{S_{1}S_{2}}(t)\\ &\quad+q_{n,1}^{U}(t)(1-\gamma_{2})p_{n}^{S_{1}I_{2}}(t),\\ p_{n}^{S_{1}R_{2}}(t+1)&=q_{n,1}^{U}(t)\gamma_{2}p_{n}^{S_{1}I_{2}}(t)+q_{n,1}^{P}(t)p_{n}^{S_{1}R_{2}}(t),\\ p_{n}^{I_{1}S_{2}}(t+1)&=[1-q_{n,1}^{U}(t)]q_{n,2}^{U}(t)p_{n}^{S_{1}S_{2}}(t)\\ &\quad+(1-\gamma_{1})q_{n,2}^{U}(t)p_{n}^{I_{1}S_{2}}(t),\\ p_{n}^{I_{1}I_{2}}(t+1)&=[1-q_{n,1}^{U}(t)][1-q_{n,2}^{U}(t)]p_{n}^{S_{1}S_{2}}(t)\\ &\quad+[1-q_{n,1}^{U}(t)](1-\gamma_{2})p_{n}^{S_{1}I_{2}}(t)\\ &\quad+(1-\gamma_{1})[1-q_{n,2}^{U}(t)]p_{n}^{I_{1}S_{2}}(t)\\ &\quad+(1-\gamma_{1})(1-\gamma_{2})p_{n}^{I_{1}I_{2}}(t),\\ p_{n}^{I_{1}R_{2}}(t+1)&=[1-q_{n,1}^{U}(t)]\gamma_{2}p_{n}^{S_{1}I_{2}}(t)\\ &\quad+[1-q_{n,1}^{P}(t)]p_{n}^{S_{1}R_{2}}(t)\\ &\quad+(1-\gamma_{1})\gamma_{2}p_{n}^{I_{1}I_{2}}(t)\\ &\quad+(1-\gamma_{1})p_{n}^{I_{1}R_{2}}(t),\\ p_{n}^{R_{1}S_{2}}(t+1)&=\gamma_{1}q_{n,2}^{U}(t)p_{n}^{I_{1}S_{2}}(t)+q_{n,2}^{P}(t)p_{n}^{R_{1}S_{2}}(t),\\ p_{n}^{R_{1}I_{2}}(t+1)&=\gamma_{1}[1-q_{n,2}^{U}(t)]p_{n}^{I_{1}S_{2}}(t)\\ &\quad+\gamma_{1}(1-\gamma_{2})p_{n}^{I_{1}I_{2}}(t)\\ &\quad+(1-q_{n,2}^{P}(t))p_{n}^{R_{1}S_{2}}(t)\\ &\quad+(1-\gamma_{2})p_{n}^{R_{1}I_{2}}(t),\\ p_{n}^{R_{1}R_{2}}(t+1)&=\gamma_{1}\gamma_{2}p_{n}^{I_{1}I_{2}}(t)+\gamma_{1}p_{n}^{I_{1}R_{2}}(t)\\ &\quad+\gamma_{2}p_{n}^{R_{1}I_{2}}(t)+p_{n}^{R_{1}R_{2}}(t).\end{split}$$

(4)

The fractions of susceptible, infected, and recovered nodes at time $t$ for disease $i$ is denoted by $\rho^{S_{i}}(t)$, $\rho^{I_{i}}(t)$, and $\rho^{R_{i}}(t)$, respectively.

	$\displaystyle\rho^{S_{i}}(t)=\frac{1}{N}\sum\limits_{n}p_{n}^{S_{i}}(t),$		(5)
	$\displaystyle\rho^{I_{i}}(t)=\frac{1}{N}\sum\limits_{n}p_{n}^{I_{i}}(t),$
	$\displaystyle\rho^{R_{i}}(t)=\frac{1}{N}\sum\limits_{n}p_{n}^{R_{i}}(t).$

Here $N$ is the number of nodes on the network. The fraction of uninfected nodes at the stationary state (when $t\rightarrow\infty$) is $\rho^{R_{i}}(0)+\rho^{S_{i}}(\infty)=v_{i}+\rho^{S_{i}}(\infty)$ for disease $i$. We perform extensive Monte Carlo simulations (MC) to validate the MMCA results obtained by Eqs. (4). 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. 1 is obtained by averaging 1000 MC simulations. The average $R^{2}$ is 0.96, indicating a high agreement between MMCA and MC simulation results.

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

$$\left.\dfrac{d\rho^{I_{i}}(t)}{dt}\right|_{t=0}<0,$$

(6)

$$\begin{split}\left.\dfrac{d\rho^{I_{i}}(t)}{dt}\right|_{t=0}=\rho^{I_{i}}(1)-\rho^{I_{i}}(0)=\frac{\sum\limits_{n}p_{n}^{I_{i}}(1)-p_{n}^{I_{i}}(0)}{N}.\end{split}$$

(7)

According to Eq. (1) and Eq. (4),

$$\begin{split}p_{n}^{I_{i}}(1)&=p_{n}^{I_{i}S_{j}}(1)+p_{n}^{I_{i}I_{j}}(1)+p_{n}^{I_{i}R_{j}}(1)\\ &=[p_{n}^{S_{i}S_{j}}(0)+p_{n}^{S_{i}I_{j}}(0)][1-q_{n,i}^{U}(0)]\\ &\quad+p_{n}^{S_{i}R_{j}}(0)[1-q_{n,i}^{P}(0)]+(1-\gamma_{i})p_{n}^{I_{i}}(0).\end{split}$$

(8)

Thus,

$$\begin{split}\left.\dfrac{d\rho^{I_{i}}(t)}{dt}\right|_{t=0}&=\frac{1}{N}\sum\limits_{n}\{[1-q_{n,i}^{U}(0)][p_{n}^{S_{i}S_{j}}(0)\\ &\quad+p_{n}^{S_{i}I_{j}}(0))]+[1-q_{n,i}^{P}(0)]p_{n}^{S_{i}R_{j}}(0)\\ &\quad-\gamma_{i}p_{n}^{I_{i}}(0)\}.\end{split}$$

(9)

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_{i}S_{j}}=(1-v_{i}-\frac{1}{N})(1-v_{j}-\frac{1}{N})$, $p_{n}^{S_{i}I_{j}}=(1-v_{i}-\frac{1}{N})\frac{1}{N}$, and $p_{n}^{S_{i}R_{j}}=(1-v_{i}-\frac{1}{N})v_{j}$. If $N$ is large enough, $p_{n}^{I_{i}}(0)=\frac{1}{N}\rightarrow 0$, consequently, from Eq. (2) and Eq. (3), we obtain

$$\begin{split}q_{n,i}^{P}(0)&\approx 1-\frac{\alpha_{i}\beta_{i}}{N}\sum\limits_{m}a_{mn}=1-\frac{\alpha_{i}\beta_{i}}{N}k_{n},\\ q_{n,i}^{U}(0)&\approx 1-\frac{\beta_{i}}{N}\sum\limits_{m}a_{mn}=1-\frac{\beta_{i}}{N}k_{n},\end{split}$$

(10)

where $k_{n}$ is the degree of node $n$. Thus,

$$\begin{split}\left.\dfrac{d\rho^{I_{i}}(t)}{dt}\right|_{t=0}&\approx\frac{\overline{k}\beta_{i}}{N}(1-v_{i})[(1-v_{j})+\alpha_{i}v_{j}]\\ &\quad-\frac{\gamma_{i}}{N},\end{split}$$

(11)

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

	$\displaystyle R_{v,1}$	$\displaystyle=R_{0,1}(1-v_{1})[1-(1-\alpha_{1})v_{2}]\leq 1,$		(12)
	$\displaystyle R_{v,2}$	$\displaystyle=R_{0,2}(1-v_{2})[1-(1-\alpha_{2})v_{1}]\leq 1,$

where $R_{0,1}=\frac{\beta_{1}\overline{k}}{\gamma_{1}}$ and $R_{0,2}=\frac{\beta_{2}\overline{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. (12) 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_{1}$ and $v_{2}$ at the beginning. Note that this corresponds to $R_{v,i}=R_{0,i}$ for $v_{1}=v_{2}=0$. According to Eqs. (12), in order to prevent the outbreaks of both diseases, the following conditions should be met

	$\displaystyle v_{1}$	$\displaystyle\geq 1-\frac{1}{R_{0,1}[1-(1-\alpha_{1})v_{2}]},$		(13)
	$\displaystyle v_{2}$	$\displaystyle\geq 1-\frac{1}{R_{0,2}[1-(1-\alpha_{2})v_{1}]},$

From Eqs. (13), 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

$$\tilde{v}_{i}=1-\frac{1}{R_{0,i}\alpha_{i}}.$$

(14)

In the following, we consider the common scenario when $R_{0,i}>1$ and $\tilde{v}_{i}\geq 0$. 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. 2. We adopt the same two-layer network setting as Fig. 1. As illustrated in Fig. 2, only points of ($v_{1},v_{2}$) in region II can prevent the outbreaks of both diseases. The intersection $(\hat{v}_{1},\hat{v}_{2})$ satisfies that

	$\displaystyle\hat{v}_{1}$	$\displaystyle=1-\frac{1}{R_{0,1}[1-(1-\alpha_{1})\hat{v}_{2}]},$		(15)
	$\displaystyle\hat{v}_{2}$	$\displaystyle=1-\frac{1}{R_{0,2}[1-(1-\alpha_{2})\hat{v}_{1}]}.$

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}}]$,

$$\begin{split}\hat{v}_{1}&=\frac{1}{2}+\frac{\alpha_{1}-\sqrt{\Delta}}{2\alpha_{1}(1-\alpha_{2})}+\frac{1-\alpha_{1}}{2R_{0,2}\alpha_{1}(1-\alpha_{2})}\\ &\quad-\frac{1}{2R_{0,1}\alpha_{1}},\\ \hat{v}_{2}&=\frac{1}{2}+\frac{\alpha_{2}-\sqrt{\Delta}}{2\alpha_{2}(1-\alpha_{1})}+\frac{1-\alpha_{2}}{2R_{0,1}\alpha_{2}(1-\alpha_{1})}\\ &\quad-\frac{1}{2R_{0,2}\alpha_{2}}.\end{split}$$

(16)

Next, we aim to identify the minimum vaccination cost that can prevent the outbreaks of both diseases. Denote $C_{1}$ and $C_{2}$ as the costs of vaccination per capita for disease 1 and disease 2, respectively. $C_{1}$ and $C_{2}$ are positive. The total vaccination cost is $V=C_{1}v_{1}N+C_{2}v_{2}N$. Then, we formulate the following optimization problem

		minimize		$\displaystyle V$		(17)
		subject to		$\displaystyle v_{1}\geq 1-\frac{1}{R_{0,1}[1-(1-\alpha_{1})v_{2}]},$
		$\displaystyle v_{2}\geq 1-\frac{1}{R_{0,2}[1-(1-\alpha_{2})v_{1}]},$
		$\displaystyle 0\leq v_{1}\leq 1,$
		$\displaystyle 0\leq v_{2}\leq 1.$

Since the slopes for $v_{1}$ and $v_{2}$ are both positive, the optimal vaccination program ($v_{1}^{*}$, $v_{2}^{*}$) is on the boundary of region II in Fig. 2.

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

$$v_{1}=1-\frac{1}{R_{0,1}[1-(1-\alpha_{1})v_{2}]}$$

and $v_{2}\in[\hat{v}_{2},1]$, thus,

$$V=C_{1}\{1-\frac{1}{R_{0,1}[1-(1-\alpha_{1})v_{2}]}\}N+C_{2}v_{2}N,$$

(18)

$$\frac{\mathrm{d}^{2}V}{\mathrm{d}v_{2}^{2}}=-\frac{2C_{1}(1-\alpha_{1})^{2}N}{R_{0,1}[1-(1-\alpha_{1})v_{2}]^{3}}<0,$$

(19)

Therefore, $V$ is a concave function of $v_{2}$. The optimal vaccination program ($v_{1}^{*}$, $v_{2}^{*}$) is either $(\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 ($v_{1}^{*}$, $v_{2}^{*}$) is either $(\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 determined by the following four nonnegative parameters:

	$\displaystyle L_{1}$	$\displaystyle=1-\hat{v}_{1},$	$\displaystyle L_{2}$	$\displaystyle=1-\hat{v}_{2},$		(20)
	$\displaystyle L_{3}$	$\displaystyle=\hat{v}_{1}-\tilde{v}_{1},$	$\displaystyle L_{4}$	$\displaystyle=\hat{v}_{2}-\tilde{v}_{2}.$		(21)

If $L_{1}L_{2}\geq L_{3}L_{4}$, we have

$$\frac{L_{4}}{L_{1}}\leq\frac{L_{2}+L_{4}}{L_{1}+L_{3}}\leq\frac{L_{2}}{L_{3}},$$

thus,

$$(v_{1}^{*},v_{2}^{*})=\left\{\begin{aligned} (1,\tilde{v}_{2})\quad&\frac{C_{1}}{C_{2}}<\frac{L_{4}}{L_{1}},\\ (\hat{v}_{1},\hat{v}_{2})\quad&\frac{L_{4}}{L_{1}}\leq\frac{C_{1}}{C_{2}}\leq\frac{L_{2}}{L_{3}},\\ (\tilde{v}_{1},1)\quad&\frac{C_{1}}{C_{2}}>\frac{L_{2}}{L_{3}}.\end{aligned}\right.$$

(22)

If $L_{1}L_{2}<L_{3}L_{4}$, we have

$$\frac{L_{4}}{L_{1}}>\frac{L_{2}+L_{4}}{L_{1}+L_{3}}>\frac{L_{2}}{L_{3}},$$

thus,

$$(v_{1}^{*},v_{2}^{*})=\left\{\begin{aligned} (1,\tilde{v}_{2})\quad&\frac{C_{1}}{C_{2}}\leq\frac{L_{2}+L_{4}}{L_{1}+L_{3}},\\ (\tilde{v}_{1},1)\quad&\text{otherwise}.\end{aligned}\right.$$

(23)

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. (22) and Eq. (23). In Fig. 3, 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. 3(a), $R_{0,1}=10$, $R_{0,2}=8$, $\alpha_{1}=0.4$, $\alpha_{2}=0.2$, so $L_{1}L_{2}=0.0607>L_{3}L_{4}=0.0221$. As a result, there are three scenarios according to the vaccine cost ratio. In Fig. 3(b), $R_{0,1}=10$, $R_{0,2}=5$, $\alpha_{1}=0.1$, $\alpha_{2}=0.2$, so $L_{1}L_{2}=0.0943<L_{3}L_{4}=0.3201$. As a result, there are two scenarios according to the vaccine cost ratio.

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. 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 Mateus et al. (2020). With limited resources, it is important to inform the model-driven vaccination programs that can contain multiple epidemics efficiently.