Uniqueness of Nash equilibrium in vaccination games

ABSTRACT

One crucial condition for the uniqueness of Nash equilibrium set in vaccination games is that the attack ratio monotonically decreases as the vaccine coverage level increasing. We consider several deterministic vaccination models in homogeneous mixing population and in heterogeneous mixing population. Based on the final size relations obtained from the deterministic epidemic models, we prove that the attack ratios can be expressed in terms of the vaccine coverage levels, and also prove that the attack ratios are decreasing functions of vaccine coverage levels. Some thresholds are presented, which depend on the vaccine efficacy. It is proved that for vaccination games in homogeneous mixing population, there is a unique Nash equilibrium for each game.

KEYWORDS:

• Epidemiology
• game theory
• uniqueness of Nash equilibrium
• vaccination game
• homogeneous mixing population
• heterogeneous mixing population
• compartmentalepidemic models
• final size relations
• population dynamics

2010 MATHEMATICS SUBJECT CLASSIFICATION:

• 92D30
• 91A80
• 92D25
• 49K15

1. Introduction

Game theory has been applied to predict human behaviour in the context of epidemiology for a long time. People choose the best strategy to maximize their own payoffs, based on the outcomes of different strategies. Before the outbreak of the epidemic, every individual can choose vaccination, or undertake the risk of infection. Ultimately, a certain percentage of people goes into the vaccinated class, while others remain in the susceptible class. Under these circumstances, the Nash equilibrium is reached; anyone adopting a different strategy will have a lower payoff.

The first paper to discuss game theory and epidemiology is [14] by Fine and Clarkson in 1986. They discussed how individuals decide whether or not to accept vaccination based on their perception of the risks involved. The conclusion of this paper is that rational informed individuals choose a lower vaccine uptake than would the community if it acted as a whole. The topic was also discussed independently by Brito et al. in [10] in 1991. Bauch and Earn described a vaccination game in their PNAS paper [3]. They characterized the Nash equilibrium of the game, proved that the Nash equilibrium is convergently stable; then they used the epidemiological model to make more precise predictions of the expected vaccine coverage level. In many publications, it has been mentioned that if individuals are driven solely by self-interest, the vaccine coverage level of the whole population will be suboptimal, and this level is lower than the herd immunity threshold value (see [3,4,15,18,19,21]). It was also mentioned that in the simple vaccination game for an endemic disease in a homogeneous mixing population, there is a unique Nash equilibrium behaviour in [18] by Reluga. One crucial condition for uniqueness is that the eventual attack ratio of an individual must decrease strictly with the vaccine coverage level [3]. The monotonicity of attack ratio is our focus in this paper. For the case of heterogeneous mixing population, in [18], Reluga considered an epidemiological game with two interacting sub-populations of equal size. He proved that if the mixing is positively assortative, there is a unique stable Nash equilibrium, while if the mixing is disassortative, this may lead to a condition of multiple Nash equilibria. It was proved that the interior Nash equilibria are always unstable. The mixing pattern of this heterogeneous mixing population is the key factor [18]. Chen also pointed out the potential for non-uniqueness in [11].

Along with the growth of interest in game theory and epidemiology, game theory has been applied to many types of infectious diseases. In [22], from the viewpoint of parents, how to make decision on Measles–Mumps–Rubella (MMR) vaccination for their children? Although scientific evidences showed that the MMR vaccine is very safe [12], there are still vaccine sceptics. There are other publications discussing the strategies against influenza [13,15,20,24], smallpox [4], chickenpox [16], measles [22], rubella [23] and other infectious diseases. Other aspects or factors of epidemiology-game theory are also discussed in many publications.

The analysis in this paper is established on the final size relations of several deterministic epidemic models. Final size relation is a very useful tool to predict how serious an epidemic can be, and can be obtained from the epidemic model. The final size relation to the simple susceptible-infected-recovered (SIR) model is $\begin{array}{rl}\ln \frac{S_0}{S_{\mathrm{\infty }}}& =\beta {\int }_0^{\mathrm{\infty }}I\left(t\right)\mathrm{d}t\\ & =\frac{\beta }{\alpha }[N-S_{\mathrm{\infty }}]\\ & ={\mathcal{R}}_0\left[1-\frac{S_{\mathrm{\infty }}}{N}\right],\end{array}$ with $S_0$ being the initial size of the susceptible class, N size of the whole population, ${\mathcal{R}}_0$ the basic reproduction number. More details can be found in papers [5,6]. This relation can be used to estimate $S_{\mathrm{\infty }}$ and the attack ratio. The same technique can be used to derive the final size relations to more complicated deterministic epidemic models. A final size relation was derived for a general class of deterministic models which includes multiple susceptible classes in [1] by Arino et al. For age-of-infection epidemic models, Brauer formulated the final size relation to this type of models in [7]. Brauer also formulated the final size relation to the epidemic model with heterogeneous mixing and treatment in [8]. For the most general age-of-infection models with heterogeneous mixing, Brauer and Watmough formulated the final size relations in [9].

The paper consists of five parts. In Section 2, we demonstrate two vaccination games in homogeneous mixing population and show why monotonicity of attack ratio is important in the analysis. Vaccination games in heterogeneous mixing population are more complicated to analyse completely. In Section 3, the focus is on the homogeneous mixing population, we discuss the relation of $S_{\mathrm{\infty }}$, $V_{\mathrm{\infty }}$ and vaccine coverage level p, and the relation between attack ratios and p, some thresholds are discussed. This part of work is able to prove the uniqueness of Nash equilibria in games which we constructed in Section 2. In Section 4, two-subgroup heterogeneous mixing population is analysed, similar analysis is conducted. We use notations $S_i\left(\mathrm{\infty }\right)$ and $V_i\left(\mathrm{\infty }\right)$ in the analysis of heterogeneous mixing population, with the subscript i indicating the subgroup in the population. In Section 5, some possible extensions are discussed. Because of the similarity of forms of final size relations, it is possible to apply the results to the general compartmental deterministic epidemic models. Finally, we propose several possible future work.

2. Vaccination games in homogeneous population

The basic method to predict the expected vaccine coverage level by game theory is coupling game theory and epidemiological models. From the analysis of game theory, an implicit form of the expected vaccine coverage level and the expected attack ratio can be obtained. Through the deterministic epidemic models, the expected attack ratio can be expressed in terms of the vaccine coverage level. Finally, an accurate estimation of vaccine coverage level can be made. We assume that the population is well-mixed/homogeneous and all individuals are provided with the same information and use this information in the same way to assess risks, this is described in [3].

2.1. Vaccination game with perfect vaccine

We formalize the vaccination game with perfect vaccine by setting the elements of the following payoff matrix:

Table

	Fully effective vaccination game
	Vaccinated individual j	Susceptible individual j
Vaccinated individual i	$-C_v$	$-C_v$
Susceptible individual i	0	$-{\pi }_p{C}_i$

This payoff matrix describes four outcomes of the game: individual i is vaccinated/not vaccinated and interacts with individual j who is vaccinated/not vaccinated. The elements in the payoff matrix are the payoffs to individual i. If individual i is vaccinated, the payoff is the cost of the vaccination $-C_v$. This is true whether this individual is matched with a vaccinated or not vaccinated individual, because the vaccine is fully effective and there is no chance for individual i to get infected. If individual i chooses not to be vaccinated, there are two possibilities. The first possibility is to match with a vaccinated individual. In this case, since individual j is protected by the vaccine and will not be infected, individual i will not be infected either, so the payoff is 0. The second possibility is to match with a not vaccinated individual. We assume for everyone who is susceptible, the probability of being infected is ${\pi }_p$, and the cost of infection is $-C_i$. We also assume that if individual j is infected, individual i will be infected through the contact. The costs of vaccination and infection were indicated as negative terms, because they decrease the payoffs. A similar game can be found in [25]. It is noted that if the cost of being vaccinated ($C_v$) is larger than the cost of being infected ($C_i$), the strategy of not being vaccinated is the dominant strategy. Since to be vaccinated is more risky than being infected.

The second step is to calculate the payoffs of individual i choosing different strategies. Assume that the final fraction of individuals who takes vaccination is p; the fraction of individuals to be vaccinated is $x_v$ and the fraction of individuals not to be vaccinated is $x_s$. $\{x_v,x_s\}=\{p,1-p\}$ is the equilibrium (or evolutionary stable strategy). We have the payoffs of two strategies: (1) $\begin{array}{rl}{f}_v& =-C_v{x}_v-C_v{x}_s=-C_v,\\ f_s& =-{\pi }_p{C}_i{x}_s,\end{array}$(1) where $f_v$ is the payoff for individuals to be vaccinated and $f_s$ is the payoff for individuals not to be vaccinated.

Then, the replicator equations to this vaccination game can be expressed as (2) $\begin{array}{rl}\dot{x_v}& =x_v(f_v-\overline{f}),\\ \dot{x_s}& =x_s(f_s-\overline{f}),\end{array}$(2) where $\overline{f}$ denotes the average payoff. By solving the replicator equations, we find that the equilibrium is $x_v=1-C_v/{\pi }_p{C}_i$, which implies that the final fraction p of population to be vaccinated is $1-C_v/{\pi }_p{C}_i$. If we define the relative cost $r=C_v/C_i$ (the relative cost reflects the relative risk of being vaccinated and being infected), we can rewrite the expected vaccine coverage level $p=1-r/{\pi }_p$.

Since ${\pi }_p$ is determined by the value of p, we rewrite the relation as (3) ${\pi }_p(1-p)=r.$(3)

2.2. Vaccination game with imperfect vaccine

To construct the vaccination game with imperfect vaccine, we introduce a factor σ, $0<\sigma <1$, in the infection rate of vaccinated members. $\sigma =0$ means that the vaccine is perfectly effective and $\sigma =1$ means that the vaccine has no effect. We set up the game by setting the entries of the payoff matrix. If individual i is vaccinated and individual j is also vaccinated, then the payoff for individual i is the cost of vaccination plus the risk of being infected, the cost of vaccination is still $-C_v$; the probability that individual j getting infected is ${\pi }_v$, since individual i is also vaccinated, the probability of being infected by individual j is reduced by σ, and the cost of infection is $-C_i$. The other three entries in the payoff matrix can be generated in the same way. It is noted that ${\pi }_p$ is the probability that individual is not vaccinated but infected. The formal partial-effective vaccination game is:

Table

	Partial-effective vaccination game
	Vaccinated individual	Susceptible individual
Vaccinated individual	$-C_v-\sigma {\pi }_v{C}_i$	$-C_v-\sigma {\pi }_p{C}_i$
Susceptible individual	$-{\pi }_v{C}_i$	$-{\pi }_p{C}_i$

The replicator equations can still be used to describe the dynamics of the game. We have the expression of the fraction p (4) $\frac{r}{1-\sigma }=p{\pi }_v+(1-p){\pi }_p,$(4) where r represents the relative cost $C_v/C_i$, ${\pi }_p$ is the attack ratio for the susceptible class and ${\pi }_v$ is the attack ratio for the vaccinated class.

2.3. Importance of property of monotonicity

In the vaccination game with perfect vaccine, we have an implicit equation (3(3) ${\pi }_p(1-p)=r.$(3) ) of the vaccine coverage level. When the vaccine coverage level p=0, the left side of the equation is ${\pi }_0$ and this value should be larger than r. When the vaccine coverage level $p\ge p_c$, $p_c$ is the herd immunity threshold value, the left side of the equation is 0, which is smaller than r. If we differentiate the function in the left side with respect to p, the derivative is ${\pi }_p^{\mathrm{\prime }}(1-p)-{\pi }_p.$ We first assume that if ${\pi }_p^{\mathrm{\prime }}=0$. In this case, ${\pi }_0={\pi }_p$ for all $0<p\le 1$. Apparently, we have ${\pi }_0>0$, which implies ${\pi }_{p_c}>0$ ($p_c$ is the herd immunity threshold). However, if $p\ge p_c$, the attack ratio has to be 0. This contradiction gives us ${\pi }_p^{\mathrm{\prime }}\ne 0$. The same method can show that it is impossible for ${\pi }_p^{\mathrm{\prime }}>0$. If ${\pi }_p^{\mathrm{\prime }}$ is negative, the derivative ${\pi }_p^{\mathrm{\prime }}(1-p)-{\pi }_p$ is negative, then Equation  (3(3) ${\pi }_p(1-p)=r.$(3) ) has a unique solution. Thus, it is important to show that the attack ratio is a decreasing function of the vaccine coverage level p.

In the vaccination game with imperfect vaccine, we have formula (4(4) $\frac{r}{1-\sigma }=p{\pi }_v+(1-p){\pi }_p,$(4) ) about the vaccine coverage level p. If the vaccine coverage level p=0, the right side is ${\pi }_0$, and this value should be larger than $r/(1-\sigma )$. If the vaccine coverage level $p\ge p_c$, the right side of the equation is smaller than $r/(1-\sigma )$. Because the left side of the equation is a constant, if the derivative of the right side is negative, we have one unique solution. We take the derivative of $p{\pi }_v+(1-p){\pi }_p$ with respect to p to get $\begin{array}{rl}& {\pi }_v+p{\pi }_v^{\mathrm{\prime }}-{\pi }_p+(1-p){\pi }_p^{\mathrm{\prime }}\\ & =({\pi }_v-{\pi }_p)+p{\pi }_v^{\mathrm{\prime }}+(1-p){\pi }_p^{\mathrm{\prime }}.\end{array}$ We have ${\pi }_v-{\pi }_p<0$, since the vaccinated individuals are protected by vaccine, the corresponding attack ratio is smaller. If both ${\pi }_v^{\mathrm{\prime }}$ and ${\pi }_p^{\mathrm{\prime }}$ are negative, then formula (4(4) $\frac{r}{1-\sigma }=p{\pi }_v+(1-p){\pi }_p,$(4) ) has a unique solution. It can be proved that ${\pi }_p^{\mathrm{\prime }}$ and ${\pi }_v^{\mathrm{\prime }}$ cannot be positive or zero.

Remark 2.1.

${\pi }_0>r$ is a necessary condition for these two vaccination games to have Nash equilibrium. Otherwise, we have two dominant games. The strategy of not to be vaccinated is dominant, because choosing the strategy of to be vaccinated is always more risky than being infected.

Remark 2.2.

The analysis of these two games implies that the herd immunity cannot be reached if all individuals are purely driven by self-interest. Since the solutions to Equations (3(3) ${\pi }_p(1-p)=r.$(3) ) and (4(4) $\frac{r}{1-\sigma }=p{\pi }_v+(1-p){\pi }_p,$(4) ) are absolutely between 0 and the herd immunity threshold $p_c$. The vaccine coverage level p we discussed in this paper is lower than $p_c$. Once we have $p\ge p_c$, the attack ratio becomes 0.

We now define the attack ratio in the population under study. This definition will be used through the paper repeatedly.

Definition 2.1.

The attack ratio in the population is defined as $\pi =\frac{N_i-N_f}{N_i}=1-\frac{N_f}{N_i},$ where $N_i$ is the size of the whole population and $N_f$ is the number of individuals who are not affected by the epidemic.

3. Mathematical analysis in homogeneous mixing population

Before the outbreak of the epidemic, people can choose to be vaccinated, or not to be vaccinated, based on the estimation of payoffs of such behaviours. We assume that the population is homogeneous mixing with size N and the expected vaccine coverage level is p, σ is the vaccine efficacy factor, with $0<\sigma <1$. $\sigma =0$ means that the vaccine is perfectly effective and $\sigma =1$ means that the vaccine has no effect. We also assume that an average member of the population makes contact sufficient to transmit infection with $\beta N$ others per unit time. Infectives leave the infective class at rate α per unit time (see [5,6] for assumptions of simple compartmental models in detail). All individuals remain in the susceptible class or go to the vaccinated class before the outbreak.

3.1. Susceptible-vaccinated-infected-recovered (SVIR) compartmental model

The SVIR model we use to describe such a condition is (5) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-\beta SI,\\ \frac{\mathrm{d}V}{\mathrm{d}t}& =-\sigma \beta VI,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\beta SI+\sigma \beta VI-\alpha I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\alpha I,\end{array}$(5) with initial values $S_0=(1-p)N$, $V_0=pN$, $I_0=0$, other parameters are defined as before. Since the population is not completely susceptible, it is accurate to use the terminology control reproduction number instead of basic reproduction number. In this case, ${\mathcal{R}}_c=\frac{\beta }{\alpha }(S_0+\sigma V_0),$ and the final size relation can be expressed by the following equations (a special case of Equation (3) in [9] by Brauer and Watmough): (6) $\begin{array}{rl}\ln \frac{S_{\mathrm{\infty }}}{S_0}& =\frac{\beta }{\alpha }[S_{\mathrm{\infty }}+V_{\mathrm{\infty }}-S_0-V_0],\\ \ln \frac{V_{\mathrm{\infty }}}{V_0}& =\frac{\sigma \beta }{\alpha }[S_{\mathrm{\infty }}+V_{\mathrm{\infty }}-S_0-V_0].\end{array}$(6) It is easy to find out that (7) $\frac{V_{\mathrm{\infty }}}{V_0}={\left(\frac{S_{\mathrm{\infty }}}{S_0}\right)}^{\sigma },$(7) which implies that for every individual, staying in group V is much safer than staying in group S.

3.2. $S_{\mathrm{\infty }}$, $V_{\mathrm{\infty }}$, attack ratios and vaccine coverage level

If we change the vaccine coverage level p, how this will affect $S_{\mathrm{\infty }}$ and $V_{\mathrm{\infty }}$? We take the derivatives of $S_{\mathrm{\infty }}$ and $V_{\mathrm{\infty }}$ with respect to p to have (8) $\begin{array}{rl}\left(\frac{1}{S_{\mathrm{\infty }}}-\frac{\beta }{\alpha }\right)\frac{\mathrm{d}{S}_{\mathrm{\infty }}}{\mathrm{d}p}& =\frac{\beta }{\alpha }\frac{\mathrm{d}{V}_{\mathrm{\infty }}}{\mathrm{d}p}-\frac{1}{1-p},\\ \left(\frac{1}{V_{\mathrm{\infty }}}-\frac{\sigma \beta }{\alpha }\right)\frac{\mathrm{d}{V}_{\mathrm{\infty }}}{\mathrm{d}p}& =\frac{\sigma \beta }{\alpha }\frac{\mathrm{d}{S}_{\mathrm{\infty }}}{\mathrm{d}p}+\frac{1}{p}.\end{array}$(8) We reorganize Equations (8(8) $\begin{array}{rl}\left(\frac{1}{S_{\mathrm{\infty }}}-\frac{\beta }{\alpha }\right)\frac{\mathrm{d}{S}_{\mathrm{\infty }}}{\mathrm{d}p}& =\frac{\beta }{\alpha }\frac{\mathrm{d}{V}_{\mathrm{\infty }}}{\mathrm{d}p}-\frac{1}{1-p},\\ \left(\frac{1}{V_{\mathrm{\infty }}}-\frac{\sigma \beta }{\alpha }\right)\frac{\mathrm{d}{V}_{\mathrm{\infty }}}{\mathrm{d}p}& =\frac{\sigma \beta }{\alpha }\frac{\mathrm{d}{S}_{\mathrm{\infty }}}{\mathrm{d}p}+\frac{1}{p}.\end{array}$(8) ) to obtasimilar equation systems (by different substitutions) (9) $\begin{array}{rl}-\frac{1}{N}\frac{\mathrm{d}{S}_{\mathrm{\infty }}}{\mathrm{d}p}& =\frac{\frac{\beta }{\alpha }\frac{\mathrm{d}{V}_{\mathrm{\infty }}}{\mathrm{d}p}(-\frac{1}{N})+\frac{1}{S_0}}{\frac{1}{S_{\mathrm{\infty }}}-\frac{\beta }{\alpha }},\\ -\frac{1}{N}\frac{1-\frac{\beta }{\alpha }{S}_{\mathrm{\infty }}-\frac{\sigma \beta }{\alpha }{V}_{\mathrm{\infty }}}{S_{\mathrm{\infty }}{V}_{\mathrm{\infty }}}\frac{\mathrm{d}{V}_{\mathrm{\infty }}}{\mathrm{d}p}& =\frac{[\frac{\sigma \beta }{\alpha }{V}_0+\frac{\beta }{\alpha }{S}_0]-\frac{S_0}{S_{\mathrm{\infty }}}}{S_0{V}_0},\end{array}$(9) and (10) $\begin{array}{rl}-\frac{1}{N}\frac{\mathrm{d}{V}_{\mathrm{\infty }}}{\mathrm{d}p}& =\frac{\frac{\sigma \beta }{\alpha }\frac{\mathrm{d}{S}_{\mathrm{\infty }}}{\mathrm{d}p}(-\frac{1}{N})-\frac{1}{V_0}}{\frac{1}{V_{\mathrm{\infty }}}-\frac{\sigma \beta }{\alpha }},\\ -\frac{1}{N}\frac{1-\frac{\sigma \beta }{\alpha }{V}_{\mathrm{\infty }}-\frac{\beta }{\alpha }{S}_{\mathrm{\infty }}}{S_{\mathrm{\infty }}}\frac{\mathrm{d}{S}_{\mathrm{\infty }}}{\mathrm{d}p}& =\frac{V_0-V_{\mathrm{\infty }}[\frac{\sigma \beta }{\alpha }{V}_0+\frac{\beta }{\alpha }{S}_0]}{S_0{V}_0}.\end{array}$(10) It is noted that in Equations (9(9) $\begin{array}{rl}-\frac{1}{N}\frac{\mathrm{d}{S}_{\mathrm{\infty }}}{\mathrm{d}p}& =\frac{\frac{\beta }{\alpha }\frac{\mathrm{d}{V}_{\mathrm{\infty }}}{\mathrm{d}p}(-\frac{1}{N})+\frac{1}{S_0}}{\frac{1}{S_{\mathrm{\infty }}}-\frac{\beta }{\alpha }},\\ -\frac{1}{N}\frac{1-\frac{\beta }{\alpha }{S}_{\mathrm{\infty }}-\frac{\sigma \beta }{\alpha }{V}_{\mathrm{\infty }}}{S_{\mathrm{\infty }}{V}_{\mathrm{\infty }}}\frac{\mathrm{d}{V}_{\mathrm{\infty }}}{\mathrm{d}p}& =\frac{[\frac{\sigma \beta }{\alpha }{V}_0+\frac{\beta }{\alpha }{S}_0]-\frac{S_0}{S_{\mathrm{\infty }}}}{S_0{V}_0},\end{array}$(9) ) and (10(10) $\begin{array}{rl}-\frac{1}{N}\frac{\mathrm{d}{V}_{\mathrm{\infty }}}{\mathrm{d}p}& =\frac{\frac{\sigma \beta }{\alpha }\frac{\mathrm{d}{S}_{\mathrm{\infty }}}{\mathrm{d}p}(-\frac{1}{N})-\frac{1}{V_0}}{\frac{1}{V_{\mathrm{\infty }}}-\frac{\sigma \beta }{\alpha }},\\ -\frac{1}{N}\frac{1-\frac{\sigma \beta }{\alpha }{V}_{\mathrm{\infty }}-\frac{\beta }{\alpha }{S}_{\mathrm{\infty }}}{S_{\mathrm{\infty }}}\frac{\mathrm{d}{S}_{\mathrm{\infty }}}{\mathrm{d}p}& =\frac{V_0-V_{\mathrm{\infty }}[\frac{\sigma \beta }{\alpha }{V}_0+\frac{\beta }{\alpha }{S}_0]}{S_0{V}_0}.\end{array}$(10) ), ${\mathcal{R}}_c=\sigma \beta /\alpha V_0+\beta /\alpha S_0$, so the key to find whether the derivatives are positive or negative is the relation of $V_{\mathrm{\infty }}/V_0$, $S_{\mathrm{\infty }}/S_0$ and $R_c$.

We first consider the condition p=0. The SVIR model is simplified to the simple SIR model, the final size relation of such model is $\ln \frac{S_0}{S_{\mathrm{\infty }}}={\mathcal{R}}_c\left[1-\frac{S_{\mathrm{\infty }}}{S_0}\right],$ where ${\mathcal{R}}_c=\beta /\alpha N$. To prove that there is a unique solution of the final size relation, we define $g\left(x\right)=\ln \frac{S_0}{x}-{\mathcal{R}}_c\left[1-\frac{x}{S_0}\right],$ we have $g\left(0_{+}\right)>0$, $g\left(S_0\right)<0$, and $g^{\mathrm{\prime }}\left(x\right)<0$ if and only if $0<x<S_0/{\mathcal{R}}_c$. The disease can only spread if ${\mathcal{R}}_c>1$. In such a case, $g\left(x\right)$ is monotone decreasing from a positive value at $x=0_{+}$ to a minimum at $x=S_0/{\mathcal{R}}_c$ and then increases to a negative value at $x=S_0$. Thus, there is a unique zero $S_{\mathrm{\infty }}$ of $g\left(x\right)$ with $S_{\mathrm{\infty }}<S_0/{\mathcal{R}}_c$, which is (11) $\frac{S_0}{S_{\mathrm{\infty }}}>{\mathcal{R}}_c,p=0.$(11)

Secondly, we consider the case that p=1. In such a case, we have a SVIR model, the final size relation of such model is $\ln \frac{V_0}{V_{\mathrm{\infty }}}={\mathcal{R}}_c\left[1-\frac{V_{\mathrm{\infty }}}{V_0}\right],$ where ${\mathcal{R}}_c=\sigma \beta /\alpha N$. If we use the same method before to prove the uniqueness of $V_{\mathrm{\infty }}$, we have (12) $\frac{V_0}{V_{\mathrm{\infty }}}>{\mathcal{R}}_c,p=1.$(12)

Now, we consider the case that 0<p<1. If we have $S_0/S_{\mathrm{\infty }}<{\mathcal{R}}_c$, the second equation in Equation (9(9) $\begin{array}{rl}-\frac{1}{N}\frac{\mathrm{d}{S}_{\mathrm{\infty }}}{\mathrm{d}p}& =\frac{\frac{\beta }{\alpha }\frac{\mathrm{d}{V}_{\mathrm{\infty }}}{\mathrm{d}p}(-\frac{1}{N})+\frac{1}{S_0}}{\frac{1}{S_{\mathrm{\infty }}}-\frac{\beta }{\alpha }},\\ -\frac{1}{N}\frac{1-\frac{\beta }{\alpha }{S}_{\mathrm{\infty }}-\frac{\sigma \beta }{\alpha }{V}_{\mathrm{\infty }}}{S_{\mathrm{\infty }}{V}_{\mathrm{\infty }}}\frac{\mathrm{d}{V}_{\mathrm{\infty }}}{\mathrm{d}p}& =\frac{[\frac{\sigma \beta }{\alpha }{V}_0+\frac{\beta }{\alpha }{S}_0]-\frac{S_0}{S_{\mathrm{\infty }}}}{S_0{V}_0},\end{array}$(9) ) yields $\mathrm{d}{V}_{\mathrm{\infty }}/\mathrm{d}p<0$, then from the first equation in Equation (9(9) $\begin{array}{rl}-\frac{1}{N}\frac{\mathrm{d}{S}_{\mathrm{\infty }}}{\mathrm{d}p}& =\frac{\frac{\beta }{\alpha }\frac{\mathrm{d}{V}_{\mathrm{\infty }}}{\mathrm{d}p}(-\frac{1}{N})+\frac{1}{S_0}}{\frac{1}{S_{\mathrm{\infty }}}-\frac{\beta }{\alpha }},\\ -\frac{1}{N}\frac{1-\frac{\beta }{\alpha }{S}_{\mathrm{\infty }}-\frac{\sigma \beta }{\alpha }{V}_{\mathrm{\infty }}}{S_{\mathrm{\infty }}{V}_{\mathrm{\infty }}}\frac{\mathrm{d}{V}_{\mathrm{\infty }}}{\mathrm{d}p}& =\frac{[\frac{\sigma \beta }{\alpha }{V}_0+\frac{\beta }{\alpha }{S}_0]-\frac{S_0}{S_{\mathrm{\infty }}}}{S_0{V}_0},\end{array}$(9) ), this leads to $\mathrm{d}{S}_{\mathrm{\infty }}/\mathrm{d}p<0$; the second equation in Equation (10(10) $\begin{array}{rl}-\frac{1}{N}\frac{\mathrm{d}{V}_{\mathrm{\infty }}}{\mathrm{d}p}& =\frac{\frac{\sigma \beta }{\alpha }\frac{\mathrm{d}{S}_{\mathrm{\infty }}}{\mathrm{d}p}(-\frac{1}{N})-\frac{1}{V_0}}{\frac{1}{V_{\mathrm{\infty }}}-\frac{\sigma \beta }{\alpha }},\\ -\frac{1}{N}\frac{1-\frac{\sigma \beta }{\alpha }{V}_{\mathrm{\infty }}-\frac{\beta }{\alpha }{S}_{\mathrm{\infty }}}{S_{\mathrm{\infty }}}\frac{\mathrm{d}{S}_{\mathrm{\infty }}}{\mathrm{d}p}& =\frac{V_0-V_{\mathrm{\infty }}[\frac{\sigma \beta }{\alpha }{V}_0+\frac{\beta }{\alpha }{S}_0]}{S_0{V}_0}.\end{array}$(10) ) yields $V_0/V_{\mathrm{\infty }}>{\mathcal{R}}_c$, finally this assumption leads to $V_0/V_{\mathrm{\infty }}>S_0/S_{\mathrm{\infty }}$. This is a contradiction to $S_0/S_{\mathrm{\infty }}>V_0/V_{\mathrm{\infty }}$. So $S_0/S_{\mathrm{\infty }}>{\mathcal{R}}_c$ must hold.

To investigate that how $S_0/S_{\mathrm{\infty }}$ and $V_0/V_{\mathrm{\infty }}$ will change if p is changed, we still use the final size relation (6(6) $\begin{array}{rl}\ln \frac{S_{\mathrm{\infty }}}{S_0}& =\frac{\beta }{\alpha }[S_{\mathrm{\infty }}+V_{\mathrm{\infty }}-S_0-V_0],\\ \ln \frac{V_{\mathrm{\infty }}}{V_0}& =\frac{\sigma \beta }{\alpha }[S_{\mathrm{\infty }}+V_{\mathrm{\infty }}-S_0-V_0].\end{array}$(6) ). For simplicity, we denote ${\pi }_1:=S_{\mathrm{\infty }}/S_0$ and ${\pi }_2:=V_{\mathrm{\infty }}/V_0$ (13) $\begin{array}{rl}\ln {\pi }_1& =\frac{\beta }{\alpha }(1-p)N\left[{\pi }_1+\frac{p}{1-p}{\pi }_2-\frac{1}{1-p}\right],\\ \ln {\pi }_2& =\frac{\sigma \beta }{\alpha }pN\left[{\pi }_2+\frac{1-p}{p}{\pi }_1-\frac{1}{p}\right].\end{array}$(13) We take the derivatives of ${\pi }_1$ and ${\pi }_2$ with respect to p to obtain (14) $\begin{array}{rl}\sigma \frac{{\pi }_{1,p}}{{\pi }_1}& =\frac{{\pi }_{2,p}}{{\pi }_2},\\ \left(1-\frac{\beta }{\alpha }{S}_{\mathrm{\infty }}-\frac{\sigma \beta }{\alpha }{V}_{\mathrm{\infty }}\right)\frac{{\pi }_{1,p}}{{\pi }_1}& =\frac{\beta }{\alpha }N({\pi }_2-{\pi }_1).\end{array}$(14) Since ${\pi }_2>{\pi }_1$, we have ${\pi }_{1,p}>0$ and ${\pi }_{2,p}>0$. From the fundamental knowledge of calculus, we know that $\frac{\mathrm{d}\left(\frac{S_0}{S_{\mathrm{\infty }}}\right)}{\mathrm{d}p}<0,\frac{\mathrm{d}\left(\frac{V_0}{V_{\mathrm{\infty }}}\right)}{\mathrm{d}p}<0,0<p<1.$ Now, we know that both $S_0/S_{\mathrm{\infty }}$ and $V_0/V_{\mathrm{\infty }}$ are monotone decreasing as p increases from 0 to 1, and for all p, we have $S_0/S_{\mathrm{\infty }}>{\mathcal{R}}_c$. If when $p=0_{+}$, $V_0/V_{\mathrm{\infty }}>{\mathcal{R}}_c$, then for all p, we have $V_0/V_{\mathrm{\infty }}>{\mathcal{R}}_c$; else if when $p=0_{+}$, $V_0/V_{\mathrm{\infty }}<{\mathcal{R}}_c$, as p increases, $V_0/V_{\mathrm{\infty }}$ will be first smaller than ${\mathcal{R}}_c$, then will be bigger than ${\mathcal{R}}_c$.

Remark 3.1.

I would like to point out, as the reviewer suggested, it is more intuitive to use the notations $S_{\mathrm{\infty }}/S_0$ and $V_{\mathrm{\infty }}/V_0$. We use the notations $S_0/S_{\mathrm{\infty }}$ and $V_0/V_{\mathrm{\infty }}$ in the analysis to simplify the calculations, because we assume that the control reproduction number ${\mathcal{R}}_c>1$.

We have two lemmas.

Lemma 3.1.

In the SVIR model (5(5) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-\beta SI,\\ \frac{\mathrm{d}V}{\mathrm{d}t}& =-\sigma \beta VI,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\beta SI+\sigma \beta VI-\alpha I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\alpha I,\end{array}$(5) ) with infection factor σ, if ${\mathcal{R}}_c>1$, we have $\beta S_0+\sigma \beta V_0>\alpha$ and $\beta S_{\mathrm{\infty }}+\sigma \beta V_{\mathrm{\infty }}<\alpha$.

Proof.

Since the assumption is that the disease can spread among people, the disease-free equilibrium is unstable, we have ${\mathcal{R}}_c=\sigma \beta /\alpha V_0+\beta /\alpha S_0>1$. And because the endemic equilibrium is asymptotically stable, so $\beta S_{\mathrm{\infty }}+\sigma \beta V_{\mathrm{\infty }}<\alpha$.

For the simple SIR model, $S_{\mathrm{\infty }}$ can be obtained from the final size relation, we define ${\sigma }_c$ as follows: (15) ${\left(\frac{S_0}{S_{\mathrm{\infty }}}\right)}^{{\sigma }_c}={\mathcal{R}}_c=\frac{\beta }{\alpha }N.$(15)

Lemma 3.2.

In the SVIR model (5(5) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-\beta SI,\\ \frac{\mathrm{d}V}{\mathrm{d}t}& =-\sigma \beta VI,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\beta SI+\sigma \beta VI-\alpha I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\alpha I,\end{array}$(5) ) with infection factor σ, there exists a threshold value ${\sigma }_c,$ which was defined in Equation (15(15) ${\left(\frac{S_0}{S_{\mathrm{\infty }}}\right)}^{{\sigma }_c}={\mathcal{R}}_c=\frac{\beta }{\alpha }N.$(15) ). And we have that if the infection factor σ is smaller than ${\sigma }_c,$ $V_0/V_{\mathrm{\infty }}<{\mathcal{R}}_c<S_0/S_{\mathrm{\infty }}$ holds when the vaccine coverage level is low, and by increasing the vaccine coverage level p, ${\mathcal{R}}_c<V_0/V_{\mathrm{\infty }}<S_0/S_{\mathrm{\infty }}$. If the infection factor σ is bigger than ${\sigma }_c,$ the inequality ${\mathcal{R}}_c<V_0/V_{\mathrm{\infty }}<S_0/S_{\mathrm{\infty }}$ holds.

The analysis establishes the result on $S_{\mathrm{\infty }}$ and $V_{\mathrm{\infty }}$.

Theorem 3.1.

In the SVIR compartmental model (5(5) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-\beta SI,\\ \frac{\mathrm{d}V}{\mathrm{d}t}& =-\sigma \beta VI,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\beta SI+\sigma \beta VI-\alpha I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\alpha I,\end{array}$(5) ), the infection factor σ will affect $S_{\mathrm{\infty }}$ and $V_{\mathrm{\infty }}$. There are two conditions. If σ is smaller than ${\sigma }_c,$ there exists a critical value $p_0$. As vaccine coverage level p increases from 0 to $p_0$, $S_{\mathrm{\infty }}$ will decrease and $V_{\mathrm{\infty }}$ will still increase; while p increases from $p_0$ to 1, $S_{\mathrm{\infty }}$ and $V_{\mathrm{\infty }}$ both increase, where $p_0$ is the critical value we described in lemma (3.2). If σ is bigger than ${\sigma }_c,$ when the vaccine coverage level p increases, $S_{\mathrm{\infty }}$ will decrease and $V_{\mathrm{\infty }}$ will increase.

Based on the definition of attack ratio we introduced before, the attack ratios ${\pi }_s$ and ${\pi }_v$ for susceptible class and vaccinated class have the following expressions: ${\pi }_s=1-\frac{S_{\mathrm{\infty }}}{S_0},{\pi }_v=1-\frac{V_{\mathrm{\infty }}}{V_0},$ since ${\pi }_{1,p}>0$ and ${\pi }_{2,p}>0$, we know that $\mathrm{d}{\pi }_s/\mathrm{d}p<0$ and $\mathrm{d}{\pi }_v/\mathrm{d}p<0$. So we reach the result on the change of attack ratios,

Theorem 3.2.

In the SVIR compartmental model, if the vaccine coverage level p increases between 0 and the herd immunity threshold, for vaccinated group, for non-vaccinated group and for the whole population, the attack ratios decrease. The infection factor σ and the attack ratios are independent.

Remark 3.2.

If we consider the vaccine coverage level p in the interval $[0,1],$ the attack ratio is a piecewise function. The function has maximum at $p=0,$ then decreases as p increasing from 0 to the herd immunity threshold $p_c$. If $p\ge p_c,$ the attack ratio is 0.

3.3. Case of perfect vaccine

In some game-epidemiology papers, people only consider the case that the vaccine is perfect, which is $\sigma =0$. In this special case, we only need to analyse the relation between attack ratio for susceptible class and the vaccine coverage level, since the vaccinated group is fully protected. We use the simple SIR model (16) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-\beta SI,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\beta SI-\alpha I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\alpha I,\end{array}$(16) with initial values $S_0=(1-p)N,$ $I_0=0$ and $R_0=0$.

The final size relation of Equation (16(16) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-\beta SI,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\beta SI-\alpha I,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\alpha I,\end{array}$(16) ) is (17) $\ln \frac{S_0}{S_{\mathrm{\infty }}}=\frac{\beta }{\alpha }{S}_0\left[1-\frac{S_{\mathrm{\infty }}}{S_0}\right],$(17) we take the derivative of $S_{\mathrm{\infty }}$ with respect to p to have (18) $\frac{\mathrm{d}{S}_{\mathrm{\infty }}}{\mathrm{d}p}=-\frac{\frac{1}{S_0}-\frac{\beta }{\alpha }}{\frac{1}{S_{\mathrm{\infty }}}-\frac{\beta }{\alpha }}>0.$(18)

We have the following results.

Corollary 3.1.

In the SIR compartmental model, when vaccine coverage level p increases from 0 to the herd immunity threshold $p_c,$ the final size $S_{\mathrm{\infty }}$ of susceptible class will increase, the attack ratio will decrease. If $p\ge p_c,$ the whole population is protected completely by the vaccine.

4. Mathematical analysis in heterogeneous mixing population

It is difficult to prove or disprove the uniqueness of Nash equilibrium in vaccination games in heterogeneous mixing population. The analysis for uniqueness consists of the properties of functions of attack ratios in subgroups and the pattern of mixing in the population. We obtain partial results on attack ratios in this section.

We consider the population is divided into subgroups with different activity levels. We will analyse how $S_i\left(\mathrm{\infty }\right)$ and $V_i\left(\mathrm{\infty }\right)$ $(i=1,2)$ in two subgroups be changed if we increase or decrease the vaccine coverage level of each subgroup. We assume that the sizes of population in the two subgroups are $N_1$ and $N_2$. Group i members make $a_i$ contacts in unit time and that the fraction of contacts made by a member of group i that is with a member of group j is $p_{ij}$, with i,j=1,2 (see [8]).

4.1. SVIR compartmental model

The two-group SVIR vaccination model is (19) $\begin{array}{rl}{S}_1^{\mathrm{\prime }}& =-\left[a_1{p}_{11}\frac{S_1{I}_1}{N_1}+a_1{p}_{12}\frac{S_1{I}_2}{N_2}\right],\\ V_1^{\mathrm{\prime }}& =-{\sigma }_1\left[a_1{p}_{11}\frac{V_1{I}_1}{N_1}+a_1{p}_{12}\frac{V_1{I}_2}{N_2}\right],\\ I_1^{\mathrm{\prime }}& =\left[a_1{p}_{11}\frac{S_1{I}_1}{N_1}+a_1{p}_{12}\frac{S_1{I}_2}{N_2}\right]+{\sigma }_1[a_1{p}_{11}\frac{V_1{I}_1}{N_1}+a_1{p}_{12}\frac{V_1{I}_2}{N_2}]-{\alpha }_1{I}_1,\\ R_1^{\mathrm{\prime }}& ={\alpha }_1{I}_1,\\ S_2^{\mathrm{\prime }}& =-\left[a_2{p}_{21}\frac{S_2{I}_1}{N_1}+a_2{p}_{22}\frac{S_2{I}_2}{N_2}\right],\\ V_2^{\mathrm{\prime }}& =-{\sigma }_2\left[a_2{p}_{21}\frac{V_2{I}_1}{N_1}+a_2{p}_{22}\frac{V_2{I}_2}{N_2}\right],\\ I_2^{\mathrm{\prime }}& =\left[a_2{p}_{21}\frac{S_2{I}_1}{N_1}+a_2{p}_{22}\frac{S_2{I}_2}{N_2}\right]+{\sigma }_2\left[a_2{p}_{21}\frac{V_2{I}_1}{N_1}+a_2{p}_{22}\frac{V_2{I}_2}{N_2}\right]-{\alpha }_2{I}_2,\\ R_2^{\mathrm{\prime }}& ={\alpha }_2{I}_2.\end{array}$(19) The mean infective period in group i is $1/{\alpha }_i$ (i=1,2). The rate of infection in group i is ${\sigma }_i$ (i=1,2) and the vaccine coverage levels in each subgroup are $v_1$ and $v_2$. The initial condition to this epidemic model is $S_1\left(0\right)=(1-v_1)N_1$, $S_2\left(0\right)=(1-v_2)N_2$, $V_1\left(0\right)=v_1{N}_1$, $V_2\left(0\right)=v_2{N}_2$, $I_1\left(0\right)=0$ and $I_2\left(0\right)=0$.

We use the next-generation matrix approach, which is described in [26], to find the control reproduction number ${\mathcal{R}}_c$. The final size relation of Equation (19(19) $\begin{array}{rl}{S}_1^{\mathrm{\prime }}& =-\left[a_1{p}_{11}\frac{S_1{I}_1}{N_1}+a_1{p}_{12}\frac{S_1{I}_2}{N_2}\right],\\ V_1^{\mathrm{\prime }}& =-{\sigma }_1\left[a_1{p}_{11}\frac{V_1{I}_1}{N_1}+a_1{p}_{12}\frac{V_1{I}_2}{N_2}\right],\\ I_1^{\mathrm{\prime }}& =\left[a_1{p}_{11}\frac{S_1{I}_1}{N_1}+a_1{p}_{12}\frac{S_1{I}_2}{N_2}\right]+{\sigma }_1[a_1{p}_{11}\frac{V_1{I}_1}{N_1}+a_1{p}_{12}\frac{V_1{I}_2}{N_2}]-{\alpha }_1{I}_1,\\ R_1^{\mathrm{\prime }}& ={\alpha }_1{I}_1,\\ S_2^{\mathrm{\prime }}& =-\left[a_2{p}_{21}\frac{S_2{I}_1}{N_1}+a_2{p}_{22}\frac{S_2{I}_2}{N_2}\right],\\ V_2^{\mathrm{\prime }}& =-{\sigma }_2\left[a_2{p}_{21}\frac{V_2{I}_1}{N_1}+a_2{p}_{22}\frac{V_2{I}_2}{N_2}\right],\\ I_2^{\mathrm{\prime }}& =\left[a_2{p}_{21}\frac{S_2{I}_1}{N_1}+a_2{p}_{22}\frac{S_2{I}_2}{N_2}\right]+{\sigma }_2\left[a_2{p}_{21}\frac{V_2{I}_1}{N_1}+a_2{p}_{22}\frac{V_2{I}_2}{N_2}\right]-{\alpha }_2{I}_2,\\ R_2^{\mathrm{\prime }}& ={\alpha }_2{I}_2.\end{array}$(19) ) is (20) $\begin{array}{rl}\ln \frac{S_1\left(0\right)}{S_1\left(\mathrm{\infty }\right)}& =\frac{a_1{p}_{11}}{{\alpha }_1{N}_1}[N_1-S_1(\mathrm{\infty })-V_1(\mathrm{\infty }\left)\right]+\frac{a_1{p}_{12}}{{\alpha }_2{N}_2}[N_2-S_2(\mathrm{\infty })-V_2(\mathrm{\infty }\left)\right],\\ \ln \frac{V_1\left(0\right)}{V_1\left(\mathrm{\infty }\right)}& ={\sigma }_1\left[\frac{a_1{p}_{11}}{{\alpha }_1{N}_1}[N_1-S_1(\mathrm{\infty })-V_1(\mathrm{\infty }\left)\right]+\frac{a_1{p}_{12}}{{\alpha }_2{N}_2}[N_2-S_2(\mathrm{\infty })-V_2(\mathrm{\infty }\left)\right]\right],\\ \ln \frac{S_2\left(0\right)}{S_2\left(\mathrm{\infty }\right)}& =\frac{a_2{p}_{21}}{{\alpha }_1{N}_1}[N_1-S_1(\mathrm{\infty })-V_1(\mathrm{\infty }\left)\right]+\frac{a_2{p}_{22}}{{\alpha }_2{N}_2}[N_2-S_2(\mathrm{\infty })-V_2(\mathrm{\infty }\left)\right],\\ \ln \frac{V_2\left(0\right)}{V_2\left(\mathrm{\infty }\right)}& ={\sigma }_2\left[\frac{a_2{p}_{21}}{{\alpha }_1{N}_1}[N_1-S_1(\mathrm{\infty })-V_1(\mathrm{\infty }\left)\right]+\frac{a_2{p}_{22}}{{\alpha }_2{N}_2}[N_2-S_2(\mathrm{\infty })-V_2(\mathrm{\infty }\left)\right]\right].\end{array}$(20) This final size relation (20(20) $\begin{array}{rl}\ln \frac{S_1\left(0\right)}{S_1\left(\mathrm{\infty }\right)}& =\frac{a_1{p}_{11}}{{\alpha }_1{N}_1}[N_1-S_1(\mathrm{\infty })-V_1(\mathrm{\infty }\left)\right]+\frac{a_1{p}_{12}}{{\alpha }_2{N}_2}[N_2-S_2(\mathrm{\infty })-V_2(\mathrm{\infty }\left)\right],\\ \ln \frac{V_1\left(0\right)}{V_1\left(\mathrm{\infty }\right)}& ={\sigma }_1\left[\frac{a_1{p}_{11}}{{\alpha }_1{N}_1}[N_1-S_1(\mathrm{\infty })-V_1(\mathrm{\infty }\left)\right]+\frac{a_1{p}_{12}}{{\alpha }_2{N}_2}[N_2-S_2(\mathrm{\infty })-V_2(\mathrm{\infty }\left)\right]\right],\\ \ln \frac{S_2\left(0\right)}{S_2\left(\mathrm{\infty }\right)}& =\frac{a_2{p}_{21}}{{\alpha }_1{N}_1}[N_1-S_1(\mathrm{\infty })-V_1(\mathrm{\infty }\left)\right]+\frac{a_2{p}_{22}}{{\alpha }_2{N}_2}[N_2-S_2(\mathrm{\infty })-V_2(\mathrm{\infty }\left)\right],\\ \ln \frac{V_2\left(0\right)}{V_2\left(\mathrm{\infty }\right)}& ={\sigma }_2\left[\frac{a_2{p}_{21}}{{\alpha }_1{N}_1}[N_1-S_1(\mathrm{\infty })-V_1(\mathrm{\infty }\left)\right]+\frac{a_2{p}_{22}}{{\alpha }_2{N}_2}[N_2-S_2(\mathrm{\infty })-V_2(\mathrm{\infty }\left)\right]\right].\end{array}$(20) ) is also a special case of Equation $\left(3\right)$ in [9] by Brauer and Watmough.

Two equations are revealed from Equation (20(20) $\begin{array}{rl}\ln \frac{S_1\left(0\right)}{S_1\left(\mathrm{\infty }\right)}& =\frac{a_1{p}_{11}}{{\alpha }_1{N}_1}[N_1-S_1(\mathrm{\infty })-V_1(\mathrm{\infty }\left)\right]+\frac{a_1{p}_{12}}{{\alpha }_2{N}_2}[N_2-S_2(\mathrm{\infty })-V_2(\mathrm{\infty }\left)\right],\\ \ln \frac{V_1\left(0\right)}{V_1\left(\mathrm{\infty }\right)}& ={\sigma }_1\left[\frac{a_1{p}_{11}}{{\alpha }_1{N}_1}[N_1-S_1(\mathrm{\infty })-V_1(\mathrm{\infty }\left)\right]+\frac{a_1{p}_{12}}{{\alpha }_2{N}_2}[N_2-S_2(\mathrm{\infty })-V_2(\mathrm{\infty }\left)\right]\right],\\ \ln \frac{S_2\left(0\right)}{S_2\left(\mathrm{\infty }\right)}& =\frac{a_2{p}_{21}}{{\alpha }_1{N}_1}[N_1-S_1(\mathrm{\infty })-V_1(\mathrm{\infty }\left)\right]+\frac{a_2{p}_{22}}{{\alpha }_2{N}_2}[N_2-S_2(\mathrm{\infty })-V_2(\mathrm{\infty }\left)\right],\\ \ln \frac{V_2\left(0\right)}{V_2\left(\mathrm{\infty }\right)}& ={\sigma }_2\left[\frac{a_2{p}_{21}}{{\alpha }_1{N}_1}[N_1-S_1(\mathrm{\infty })-V_1(\mathrm{\infty }\left)\right]+\frac{a_2{p}_{22}}{{\alpha }_2{N}_2}[N_2-S_2(\mathrm{\infty })-V_2(\mathrm{\infty }\left)\right]\right].\end{array}$(20) ) (21) $\begin{array}{rl}{\left(\frac{S_1\left(0\right)}{S_1\left(\mathrm{\infty }\right)}\right)}^{{\sigma }_1}& =\frac{V_1\left(0\right)}{V_1\left(\mathrm{\infty }\right)},\\ {\left(\frac{S_2\left(0\right)}{S_2\left(\mathrm{\infty }\right)}\right)}^{{\sigma }_2}& =\frac{V_2\left(0\right)}{V_2\left(\mathrm{\infty }\right)}.\end{array}$(21)

4.2. $S_{\mathrm{\infty }}$, $V_{\mathrm{\infty }}$, attack ratios and vaccine coverage level

For simplicity of further calculation, we denote $c_1:=\frac{a_1{p}_{11}}{{\alpha }_1{N}_1},c_2:=\frac{a_1{p}_{12}}{{\alpha }_2{N}_2},c_3:=\frac{a_2{p}_{21}}{{\alpha }_1{N}_1},c_4:=\frac{a_2{p}_{22}}{{\alpha }_2{N}_2}.$ We take the derivatives of $S_1\left(\mathrm{\infty }\right)$, $V_1\left(\mathrm{\infty }\right)$, $S_2\left(\mathrm{\infty }\right)$ and $V_2\left(\mathrm{\infty }\right)$ with respect to $v_1$, to find out if we only adjust vaccine coverage level in one subgroup, how this will affect two subgroups. (22) $\begin{array}{rl}-\frac{1}{1-v_1}-\frac{S_{1,v_1}\left(\mathrm{\infty }\right)}{S_1\left(\mathrm{\infty }\right)}& =c_1[-S_{1,v_1}(\mathrm{\infty })-V_{1,v_1}(\mathrm{\infty }\left)\right]+c_2[-S_{2,v_1}(\mathrm{\infty })-V_{2,v_1}(\mathrm{\infty }\left)\right],\\ \frac{1}{v_1}-\frac{V_{1,v_1}\left(\mathrm{\infty }\right)}{V_1\left(\mathrm{\infty }\right)}& ={\sigma }_1{c}_1[-S_{1,v_1}(\mathrm{\infty })-V_{1,v_1}(\mathrm{\infty }\left)\right]+{\sigma }_1{c}_2[-S_{2,v_1}(\mathrm{\infty })-V_{2,v_1}(\mathrm{\infty }\left)\right],\\ -\frac{S_{2,v_1}\left(\mathrm{\infty }\right)}{S_2\left(\mathrm{\infty }\right)}& =c_3[-S_{1,v_1}(\mathrm{\infty })-V_{1,v_1}(\mathrm{\infty }\left)\right]+c_2[-S_{2,v_1}(\mathrm{\infty })-V_{2,v_1}(\mathrm{\infty }\left)\right],\\ -\frac{V_{2,v_1}\left(\mathrm{\infty }\right)}{V_2\left(\mathrm{\infty }\right)}& ={\sigma }_2{c}_3[-S_{1,v_1}(\mathrm{\infty })-V_{1,v_1}(\mathrm{\infty }\left)\right]+{\sigma }_2{c}_4[-S_{2,v_1}(\mathrm{\infty })-V_{2,v_1}(\mathrm{\infty }\left)\right],\end{array}$(22) with $S_{i,v_1}\left(\mathrm{\infty }\right):=\frac{\mathrm{d}\left(S_i\right(\mathrm{\infty }\left)\right)}{\mathrm{d}{v}_1},V_{i,v_1}\left(\mathrm{\infty }\right):=\frac{\mathrm{d}\left(V_i\right(\mathrm{\infty }\left)\right)}{\mathrm{d}{v}_1}(i=1,2).$ From Equation (22(22) $\begin{array}{rl}-\frac{1}{1-v_1}-\frac{S_{1,v_1}\left(\mathrm{\infty }\right)}{S_1\left(\mathrm{\infty }\right)}& =c_1[-S_{1,v_1}(\mathrm{\infty })-V_{1,v_1}(\mathrm{\infty }\left)\right]+c_2[-S_{2,v_1}(\mathrm{\infty })-V_{2,v_1}(\mathrm{\infty }\left)\right],\\ \frac{1}{v_1}-\frac{V_{1,v_1}\left(\mathrm{\infty }\right)}{V_1\left(\mathrm{\infty }\right)}& ={\sigma }_1{c}_1[-S_{1,v_1}(\mathrm{\infty })-V_{1,v_1}(\mathrm{\infty }\left)\right]+{\sigma }_1{c}_2[-S_{2,v_1}(\mathrm{\infty })-V_{2,v_1}(\mathrm{\infty }\left)\right],\\ -\frac{S_{2,v_1}\left(\mathrm{\infty }\right)}{S_2\left(\mathrm{\infty }\right)}& =c_3[-S_{1,v_1}(\mathrm{\infty })-V_{1,v_1}(\mathrm{\infty }\left)\right]+c_2[-S_{2,v_1}(\mathrm{\infty })-V_{2,v_1}(\mathrm{\infty }\left)\right],\\ -\frac{V_{2,v_1}\left(\mathrm{\infty }\right)}{V_2\left(\mathrm{\infty }\right)}& ={\sigma }_2{c}_3[-S_{1,v_1}(\mathrm{\infty })-V_{1,v_1}(\mathrm{\infty }\left)\right]+{\sigma }_2{c}_4[-S_{2,v_1}(\mathrm{\infty })-V_{2,v_1}(\mathrm{\infty }\left)\right],\end{array}$(22) ), we have (23) $\begin{array}{rl}\frac{V_{1,v_1}\left(\mathrm{\infty }\right)}{V_1\left(\mathrm{\infty }\right)}& ={\sigma }_1\frac{S_{1,v_1}\left(\mathrm{\infty }\right)}{S_1\left(\mathrm{\infty }\right)}+\frac{{\sigma }_1}{1-v_1}+\frac{1}{v_1},\\ \frac{V_{2,v_1}\left(\mathrm{\infty }\right)}{V_2\left(\mathrm{\infty }\right)}& ={\sigma }_2\frac{S_{2,v_1}\left(\mathrm{\infty }\right)}{S_2\left(\mathrm{\infty }\right)}.\end{array}$(23) The above relation implies that $V_{2,v_1}\left(\mathrm{\infty }\right)$ and $S_{2,v_1}\left(\mathrm{\infty }\right)$ are both positive or negative, which indicates the change of vaccine coverage level in subgroup 1 having the same impact to classes $V_2$ and $S_2$ in subgroup 2.

We state the following lemma, and use the lemma to do further analysis.

Lemma 4.1.

In the SVIR heterogeneous model (19(19) $\begin{array}{rl}{S}_1^{\mathrm{\prime }}& =-\left[a_1{p}_{11}\frac{S_1{I}_1}{N_1}+a_1{p}_{12}\frac{S_1{I}_2}{N_2}\right],\\ V_1^{\mathrm{\prime }}& =-{\sigma }_1\left[a_1{p}_{11}\frac{V_1{I}_1}{N_1}+a_1{p}_{12}\frac{V_1{I}_2}{N_2}\right],\\ I_1^{\mathrm{\prime }}& =\left[a_1{p}_{11}\frac{S_1{I}_1}{N_1}+a_1{p}_{12}\frac{S_1{I}_2}{N_2}\right]+{\sigma }_1[a_1{p}_{11}\frac{V_1{I}_1}{N_1}+a_1{p}_{12}\frac{V_1{I}_2}{N_2}]-{\alpha }_1{I}_1,\\ R_1^{\mathrm{\prime }}& ={\alpha }_1{I}_1,\\ S_2^{\mathrm{\prime }}& =-\left[a_2{p}_{21}\frac{S_2{I}_1}{N_1}+a_2{p}_{22}\frac{S_2{I}_2}{N_2}\right],\\ V_2^{\mathrm{\prime }}& =-{\sigma }_2\left[a_2{p}_{21}\frac{V_2{I}_1}{N_1}+a_2{p}_{22}\frac{V_2{I}_2}{N_2}\right],\\ I_2^{\mathrm{\prime }}& =\left[a_2{p}_{21}\frac{S_2{I}_1}{N_1}+a_2{p}_{22}\frac{S_2{I}_2}{N_2}\right]+{\sigma }_2\left[a_2{p}_{21}\frac{V_2{I}_1}{N_1}+a_2{p}_{22}\frac{V_2{I}_2}{N_2}\right]-{\alpha }_2{I}_2,\\ R_2^{\mathrm{\prime }}& ={\alpha }_2{I}_2.\end{array}$(19) ), if the epidemic can spread among the population initially, we have (24) $\begin{array}{rl}{c}_1{S}_1\left(0\right)+{\sigma }_1{c}_1{V}_1\left(0\right)& >1,\\ c_4{S}_2\left(0\right)+{\sigma }_2{c}_4{V}_2\left(0\right)& >1,\end{array}$(24) and (25) $\begin{array}{rl}{c}_1{S}_1\left(\mathrm{\infty }\right)+{\sigma }_1{c}_1{V}_1\left(\mathrm{\infty }\right)& <1,\\ c_4{S}_2\left(\mathrm{\infty }\right)+{\sigma }_2{c}_4{V}_2\left(\mathrm{\infty }\right)& <1,\end{array}$(25) $c_1,$ $c_2,$ $c_3$ and $c_4$ are defined before.

Proof.

We first linearize the system at the disease-free equilibrium to obtain $\left[\begin{array}{c}{I}_1^{\mathrm{\prime }}\\ I_2^{\mathrm{\prime }}\end{array}\right]=\left[\begin{array}{cc}{\alpha }_1\left(c_1{S}_1\right(0)+{\sigma }_1{c}_1{V}_1(0)-1)& {\alpha }_2\left(c_2{S}_1\right(0)+{\sigma }_1{c}_2{V}_1(0\left)\right)\\ {\alpha }_1\left(c_3{S}_2\right(0)+{\sigma }_2{c}_3{V}_2(0\left)\right)& {\alpha }_2\left(c_4{S}_2\right(0)+{\sigma }_2{c}_4{V}_2(0)-1)\end{array}\right]\left[\begin{array}{c}{I}_1\\ I_2\end{array}\right].$ This equilibrium is unstable, the $2\times 2$ matrix has two positive eigenvalues because the disease can spread in two subgroups. $\begin{array}{rl}& {\alpha }_1\left(c_1{S}_1\right(0)+{\sigma }_1{c}_1{V}_1(0)-1)+{\alpha }_2\left(c_4{S}_2\right(0)+{\sigma }_2{c}_4{V}_2(0)-1)>0,\\ & \left(c_1{S}_1\right(0)+{\sigma }_1{c}_1{V}_1(0)-1)\left(c_4{S}_2\right(0)+{\sigma }_2{c}_4{V}_2(0)-1)\\ & >\left(c_2{S}_1\right(0)+{\sigma }_1{c}_2{V}_1(0\left)\right)\left(c_3{S}_2\right(0)+{\sigma }_2{c}_3{V}_2(0\left)\right),\end{array}$ which lead to $c_1{S}_1\left(0\right)+{\sigma }_1{c}_1{V}_1\left(0\right)>1,$ and $c_4{S}_2\left(0\right)+{\sigma }_2{c}_4{V}_2\left(0\right)>1.$

We then linearize the system at the endemic equilibrium to obtain $\left[\begin{array}{c}{I}_1^{\mathrm{\prime }}\\ I_2^{\mathrm{\prime }}\end{array}\right]=\left[\begin{array}{cc}{\alpha }_1\left(c_1{S}_1\right(\mathrm{\infty })+{\sigma }_1{c}_1{V}_1(\mathrm{\infty })-1)& {\alpha }_2\left(c_2{S}_1\right(\mathrm{\infty })+{\sigma }_1{c}_2{V}_1(\mathrm{\infty }\left)\right)\\ {\alpha }_1\left(c_3{S}_2\right(\mathrm{\infty })+{\sigma }_2{c}_3{V}_2(\mathrm{\infty }\left)\right)& {\alpha }_2\left(c_4{S}_2\right(\mathrm{\infty })+{\sigma }_2{c}_4{V}_2(\mathrm{\infty })-1)\end{array}\right]\left[\begin{array}{c}{I}_1\\ I_2\end{array}\right].$ The endemic equilibrium is asymptotically stable, so this $2\times 2$ matrix has two negative eigenvalues $\begin{array}{rl}& {\alpha }_1\left(c_1{S}_1\right(\mathrm{\infty })+{\sigma }_1{c}_1{V}_1(\mathrm{\infty })-1)+{\alpha }_2\left(c_4{S}_2\right(\mathrm{\infty })+{\sigma }_2{c}_4{V}_2(\mathrm{\infty })-1)<0,\\ & \left(c_1{S}_1\right(\mathrm{\infty })+{\sigma }_1{c}_1{V}_1(\mathrm{\infty })-1)\left(c_4{S}_2\right(\mathrm{\infty })+{\sigma }_2{c}_4{V}_2(\mathrm{\infty })-1)\\ & >c_2{c}_3\left(S_1\right(\mathrm{\infty })+{\sigma }_1{V}_1(\mathrm{\infty }\left)\right)\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right),\end{array}$ which lead to $c_1{S}_1\left(\mathrm{\infty }\right)+{\sigma }_1{c}_1{V}_1\left(\mathrm{\infty }\right)<1,$ and $c_4{S}_2\left(\mathrm{\infty }\right)+{\sigma }_2{c}_4{V}_2\left(\mathrm{\infty }\right)<1.$

We now substitute expressions in Equation (23(23) $\begin{array}{rl}\frac{V_{1,v_1}\left(\mathrm{\infty }\right)}{V_1\left(\mathrm{\infty }\right)}& ={\sigma }_1\frac{S_{1,v_1}\left(\mathrm{\infty }\right)}{S_1\left(\mathrm{\infty }\right)}+\frac{{\sigma }_1}{1-v_1}+\frac{1}{v_1},\\ \frac{V_{2,v_1}\left(\mathrm{\infty }\right)}{V_2\left(\mathrm{\infty }\right)}& ={\sigma }_2\frac{S_{2,v_1}\left(\mathrm{\infty }\right)}{S_2\left(\mathrm{\infty }\right)}.\end{array}$(23) ) into the first and the third equation in Equation (22(22) $\begin{array}{rl}-\frac{1}{1-v_1}-\frac{S_{1,v_1}\left(\mathrm{\infty }\right)}{S_1\left(\mathrm{\infty }\right)}& =c_1[-S_{1,v_1}(\mathrm{\infty })-V_{1,v_1}(\mathrm{\infty }\left)\right]+c_2[-S_{2,v_1}(\mathrm{\infty })-V_{2,v_1}(\mathrm{\infty }\left)\right],\\ \frac{1}{v_1}-\frac{V_{1,v_1}\left(\mathrm{\infty }\right)}{V_1\left(\mathrm{\infty }\right)}& ={\sigma }_1{c}_1[-S_{1,v_1}(\mathrm{\infty })-V_{1,v_1}(\mathrm{\infty }\left)\right]+{\sigma }_1{c}_2[-S_{2,v_1}(\mathrm{\infty })-V_{2,v_1}(\mathrm{\infty }\left)\right],\\ -\frac{S_{2,v_1}\left(\mathrm{\infty }\right)}{S_2\left(\mathrm{\infty }\right)}& =c_3[-S_{1,v_1}(\mathrm{\infty })-V_{1,v_1}(\mathrm{\infty }\left)\right]+c_2[-S_{2,v_1}(\mathrm{\infty })-V_{2,v_1}(\mathrm{\infty }\left)\right],\\ -\frac{V_{2,v_1}\left(\mathrm{\infty }\right)}{V_2\left(\mathrm{\infty }\right)}& ={\sigma }_2{c}_3[-S_{1,v_1}(\mathrm{\infty })-V_{1,v_1}(\mathrm{\infty }\left)\right]+{\sigma }_2{c}_4[-S_{2,v_1}(\mathrm{\infty })-V_{2,v_1}(\mathrm{\infty }\left)\right],\end{array}$(22) ) (26) $\begin{array}{rl}& \left(\frac{1}{S_1\left(\mathrm{\infty }\right)}-c_1(1+\frac{{\sigma }_1{V}_1\left(\mathrm{\infty }\right)}{S_1\left(\mathrm{\infty }\right)})\right)S_{1,v_1}\left(\mathrm{\infty }\right)-c_2\left(1+\frac{{\sigma }_2{V}_2\left(\mathrm{\infty }\right)}{S_2\left(\mathrm{\infty }\right)}\right)S_{2,v_1}\left(\mathrm{\infty }\right)\\ & =c_1\left(\frac{{\sigma }_1{V}_1\left(\mathrm{\infty }\right)}{1-v_1}+\frac{V_1\left(\mathrm{\infty }\right)}{v_1}\right)-\frac{1}{1-v_1},\\ & \left(\frac{1}{S_2\left(\mathrm{\infty }\right)}-c_4(1+\frac{{\sigma }_2{V}_2\left(\mathrm{\infty }\right)}{S_2\left(\mathrm{\infty }\right)})\right)S_{1,v_1}\left(\mathrm{\infty }\right)-c_3\left(1+\frac{{\sigma }_1{V}_1\left(\mathrm{\infty }\right)}{S_1\left(\mathrm{\infty }\right)}\right)S_{1,v_1}\left(\mathrm{\infty }\right)\\ & =c_3\left(\frac{{\sigma }_1{V}_1\left(\mathrm{\infty }\right)}{1-v_1}+\frac{V_1\left(\mathrm{\infty }\right)}{v_1}\right).\end{array}$(26) We extract the expression of $\mathrm{\partial }{S}_2/\mathrm{\partial }{v}_1$ (27) $\begin{array}{rl}& \left[\left(\frac{1}{S_2\left(\mathrm{\infty }\right)}-c_4\left(1+{\sigma }_2\frac{V_2\left(\mathrm{\infty }\right)}{S_2\left(\mathrm{\infty }\right)}\right)\right)\left(\frac{1}{S_1\left(\mathrm{\infty }\right)}-c_1\left(1+{\sigma }_2\frac{V_1\left(\mathrm{\infty }\right)}{S_1\left(\mathrm{\infty }\right)}\right)\right)\right.\\ & \left.-c_2{c}_3(1+{\sigma }_1\frac{V_1\left(\mathrm{\infty }\right)}{S_1\left(\mathrm{\infty }\right)})\left(1+{\sigma }_1\frac{V_2\left(\mathrm{\infty }\right)}{S_2\left(\mathrm{\infty }\right)}\right)\right]\frac{\mathrm{\partial }{S}_2}{\mathrm{\partial }{v}_1}=c_3\left(\frac{V_1\left(\mathrm{\infty }\right)}{v_1{S}_1\left(\mathrm{\infty }\right)}-\frac{1}{1-v_1}\right).\end{array}$(27) Based on Equation (21(21) $\begin{array}{rl}{\left(\frac{S_1\left(0\right)}{S_1\left(\mathrm{\infty }\right)}\right)}^{{\sigma }_1}& =\frac{V_1\left(0\right)}{V_1\left(\mathrm{\infty }\right)},\\ {\left(\frac{S_2\left(0\right)}{S_2\left(\mathrm{\infty }\right)}\right)}^{{\sigma }_2}& =\frac{V_2\left(0\right)}{V_2\left(\mathrm{\infty }\right)}.\end{array}$(21) ), we know that the right-hand side of Equation (27(27) $\begin{array}{rl}& \left[\left(\frac{1}{S_2\left(\mathrm{\infty }\right)}-c_4\left(1+{\sigma }_2\frac{V_2\left(\mathrm{\infty }\right)}{S_2\left(\mathrm{\infty }\right)}\right)\right)\left(\frac{1}{S_1\left(\mathrm{\infty }\right)}-c_1\left(1+{\sigma }_2\frac{V_1\left(\mathrm{\infty }\right)}{S_1\left(\mathrm{\infty }\right)}\right)\right)\right.\\ & \left.-c_2{c}_3(1+{\sigma }_1\frac{V_1\left(\mathrm{\infty }\right)}{S_1\left(\mathrm{\infty }\right)})\left(1+{\sigma }_1\frac{V_2\left(\mathrm{\infty }\right)}{S_2\left(\mathrm{\infty }\right)}\right)\right]\frac{\mathrm{\partial }{S}_2}{\mathrm{\partial }{v}_1}=c_3\left(\frac{V_1\left(\mathrm{\infty }\right)}{v_1{S}_1\left(\mathrm{\infty }\right)}-\frac{1}{1-v_1}\right).\end{array}$(27) ) is positive. And from the second part of the proof of previous lemma, we know that the left-hand side of Equation (27(27) $\begin{array}{rl}& \left[\left(\frac{1}{S_2\left(\mathrm{\infty }\right)}-c_4\left(1+{\sigma }_2\frac{V_2\left(\mathrm{\infty }\right)}{S_2\left(\mathrm{\infty }\right)}\right)\right)\left(\frac{1}{S_1\left(\mathrm{\infty }\right)}-c_1\left(1+{\sigma }_2\frac{V_1\left(\mathrm{\infty }\right)}{S_1\left(\mathrm{\infty }\right)}\right)\right)\right.\\ & \left.-c_2{c}_3(1+{\sigma }_1\frac{V_1\left(\mathrm{\infty }\right)}{S_1\left(\mathrm{\infty }\right)})\left(1+{\sigma }_1\frac{V_2\left(\mathrm{\infty }\right)}{S_2\left(\mathrm{\infty }\right)}\right)\right]\frac{\mathrm{\partial }{S}_2}{\mathrm{\partial }{v}_1}=c_3\left(\frac{V_1\left(\mathrm{\infty }\right)}{v_1{S}_1\left(\mathrm{\infty }\right)}-\frac{1}{1-v_1}\right).\end{array}$(27) ) is also positive, which gives that (28) $\frac{\mathrm{\partial }{S}_2}{\mathrm{\partial }{v}_1}>0.$(28) The above analysis implies that if the vaccine coverage level in subgroup 1 is increased, for subgroup 2, both $S_2\left(\mathrm{\infty }\right)$ and $V_2\left(\mathrm{\infty }\right)$ will increase. We now investigate how $S_1\left(\mathrm{\infty }\right)$ and $V_1\left(\mathrm{\infty }\right)$ will be affected, and how the change will relate to the vaccine efficacy. Through some simple calculations, we have the left-hand side of the expression of $\mathrm{\partial }{S}_1/\mathrm{\partial }{v}_1$ be (29) $\begin{array}{rl}& \frac{1}{S_1\left(\mathrm{\infty }\right)S_2\left(\mathrm{\infty }\right)}[1-c_4(S_2\left(\mathrm{\infty }\right)+{\sigma }_2{V}_2\left(\mathrm{\infty }\right))-c_1(S_1\left(\mathrm{\infty }\right)+{\sigma }_1{V}_1\left(\mathrm{\infty }\right))\\ & +(c_1{c}_4-c_2{c}_3)\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)\left(S_1\right(\mathrm{\infty })+{\sigma }_1{V}_1(\mathrm{\infty }\left)\right)].\end{array}$(29) The right-hand side of the expression of $\mathrm{\partial }{S}_1\left(\mathrm{\infty }\right)/\mathrm{\partial }{v}_1$ is (30) $\begin{array}{rl}& \frac{N_1}{S_2\left(\mathrm{\infty }\right)}\left[(c_2{c}_3-c_1{c}_4)\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)\left(\frac{{\sigma }_1{V}_1\left(\mathrm{\infty }\right)}{S_1\left(0\right)}+\frac{V_1\left(\mathrm{\infty }\right)}{V_1\left(0\right)}\right)\right.\\ & \left.+c_1\left(\frac{{\sigma }_1{V}_1\left(\mathrm{\infty }\right)}{S_1\left(0\right)}+\frac{V_1\left(\mathrm{\infty }\right)}{V_1\left(0\right)}\right)+\frac{c_4}{S_1\left(0\right)}\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)-\frac{1}{S_1\left(0\right)}\right].\end{array}$(30) By observing the left-hand sides of Equations (27(27) $\begin{array}{rl}& \left[\left(\frac{1}{S_2\left(\mathrm{\infty }\right)}-c_4\left(1+{\sigma }_2\frac{V_2\left(\mathrm{\infty }\right)}{S_2\left(\mathrm{\infty }\right)}\right)\right)\left(\frac{1}{S_1\left(\mathrm{\infty }\right)}-c_1\left(1+{\sigma }_2\frac{V_1\left(\mathrm{\infty }\right)}{S_1\left(\mathrm{\infty }\right)}\right)\right)\right.\\ & \left.-c_2{c}_3(1+{\sigma }_1\frac{V_1\left(\mathrm{\infty }\right)}{S_1\left(\mathrm{\infty }\right)})\left(1+{\sigma }_1\frac{V_2\left(\mathrm{\infty }\right)}{S_2\left(\mathrm{\infty }\right)}\right)\right]\frac{\mathrm{\partial }{S}_2}{\mathrm{\partial }{v}_1}=c_3\left(\frac{V_1\left(\mathrm{\infty }\right)}{v_1{S}_1\left(\mathrm{\infty }\right)}-\frac{1}{1-v_1}\right).\end{array}$(27) ) and (29(29) $\begin{array}{rl}& \frac{1}{S_1\left(\mathrm{\infty }\right)S_2\left(\mathrm{\infty }\right)}[1-c_4(S_2\left(\mathrm{\infty }\right)+{\sigma }_2{V}_2\left(\mathrm{\infty }\right))-c_1(S_1\left(\mathrm{\infty }\right)+{\sigma }_1{V}_1\left(\mathrm{\infty }\right))\\ & +(c_1{c}_4-c_2{c}_3)\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)\left(S_1\right(\mathrm{\infty })+{\sigma }_1{V}_1(\mathrm{\infty }\left)\right)].\end{array}$(29) ), they have the same form, so we only need to analyse the positive or negative of Equation (30(30) $\begin{array}{rl}& \frac{N_1}{S_2\left(\mathrm{\infty }\right)}\left[(c_2{c}_3-c_1{c}_4)\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)\left(\frac{{\sigma }_1{V}_1\left(\mathrm{\infty }\right)}{S_1\left(0\right)}+\frac{V_1\left(\mathrm{\infty }\right)}{V_1\left(0\right)}\right)\right.\\ & \left.+c_1\left(\frac{{\sigma }_1{V}_1\left(\mathrm{\infty }\right)}{S_1\left(0\right)}+\frac{V_1\left(\mathrm{\infty }\right)}{V_1\left(0\right)}\right)+\frac{c_4}{S_1\left(0\right)}\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)-\frac{1}{S_1\left(0\right)}\right].\end{array}$(30) ), to rewrite Equation (30(30) $\begin{array}{rl}& \frac{N_1}{S_2\left(\mathrm{\infty }\right)}\left[(c_2{c}_3-c_1{c}_4)\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)\left(\frac{{\sigma }_1{V}_1\left(\mathrm{\infty }\right)}{S_1\left(0\right)}+\frac{V_1\left(\mathrm{\infty }\right)}{V_1\left(0\right)}\right)\right.\\ & \left.+c_1\left(\frac{{\sigma }_1{V}_1\left(\mathrm{\infty }\right)}{S_1\left(0\right)}+\frac{V_1\left(\mathrm{\infty }\right)}{V_1\left(0\right)}\right)+\frac{c_4}{S_1\left(0\right)}\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)-\frac{1}{S_1\left(0\right)}\right].\end{array}$(30) ) as (31) $\begin{array}{rl}& \frac{N_1}{S_2\left(\mathrm{\infty }\right)}\left[(c_2{c}_3-c_1{c}_4)\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)\left(\frac{{\sigma }_1{V}_1\left(\mathrm{\infty }\right)}{S_1\left(0\right)}+\frac{V_1\left(\mathrm{\infty }\right)}{V_1\left(0\right)}\right)\right.\\ & \left.+c_1\left(\frac{{\sigma }_1{V}_1\left(\mathrm{\infty }\right)}{S_1\left(0\right)}+\frac{V_1\left(\mathrm{\infty }\right)}{V_1\left(0\right)}\right)+\frac{c_4}{S_1\left(0\right)}\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)-\frac{1}{S_1\left(0\right)}\right]\\ & =\frac{N_1}{S_2\left(\mathrm{\infty }\right)S_1\left(0\right)}\left[(c_2{c}_3-c_1{c}_4)\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)\left(S_1\right(\mathrm{\infty })+{\sigma }_1{V}_1(\mathrm{\infty }\left)\right)\right.\\ & \left.+c_1({\sigma }_1{V}_1\left(\mathrm{\infty }\right)+S_1\left(\mathrm{\infty }\right))+c_4\left({\sigma }_2{V}_2\right(\mathrm{\infty })+S_2(\mathrm{\infty }\left)\right)-1\right]\\ & +\frac{N_1}{S_2\left(\mathrm{\infty }\right)S_1\left(0\right)}\left[\right(c_2{c}_3-c_1{c}_4\left)\right(S_1\left(\mathrm{\infty }\right)+{\sigma }_1{V}_1\left(\mathrm{\infty }\right))+c_1]\left(\frac{S_1\left(0\right)V_1\left(\mathrm{\infty }\right)}{V_1\left(0\right)}-S_1\left(\mathrm{\infty }\right)\right).\end{array}$(31) If ${\sigma }_1=1$, then $\frac{S_1\left(0\right)V_1\left(\mathrm{\infty }\right)}{V_1\left(0\right)}=S_1\left(\mathrm{\infty }\right),$ since the vaccine has no effect. In such a case, we have (32) $\frac{\mathrm{\partial }{S}_1\left(\mathrm{\infty }\right)}{\mathrm{\partial }{v}_1}=-\frac{N_1{S}_1\left(\mathrm{\infty }\right)}{S_1\left(0\right)}<0.$(32) Due to continuity, we have that when the vaccine is low-effective, ${\sigma }_1$ is close to 1, then if the vaccine coverage level in subgroup 1 is increased, $S_1\left(\mathrm{\infty }\right)$ will decrease. If we substitute Equation (32(32) $\frac{\mathrm{\partial }{S}_1\left(\mathrm{\infty }\right)}{\mathrm{\partial }{v}_1}=-\frac{N_1{S}_1\left(\mathrm{\infty }\right)}{S_1\left(0\right)}<0.$(32) ) into the first equation of Equation (23(23) $\begin{array}{rl}\frac{V_{1,v_1}\left(\mathrm{\infty }\right)}{V_1\left(\mathrm{\infty }\right)}& ={\sigma }_1\frac{S_{1,v_1}\left(\mathrm{\infty }\right)}{S_1\left(\mathrm{\infty }\right)}+\frac{{\sigma }_1}{1-v_1}+\frac{1}{v_1},\\ \frac{V_{2,v_1}\left(\mathrm{\infty }\right)}{V_2\left(\mathrm{\infty }\right)}& ={\sigma }_2\frac{S_{2,v_1}\left(\mathrm{\infty }\right)}{S_2\left(\mathrm{\infty }\right)}.\end{array}$(23) ), we have $\begin{array}{rl}\frac{1}{V_1\left(\mathrm{\infty }\right)}\frac{\mathrm{\partial }{V}_1\left(\mathrm{\infty }\right)}{\mathrm{\partial }{v}_1}& ={\sigma }_1\frac{1}{S_1\left(\mathrm{\infty }\right)}\frac{\mathrm{\partial }{S}_1\left(\mathrm{\infty }\right)}{\mathrm{\partial }{v}_1}+\frac{{\sigma }_1}{1-v_1}+\frac{1}{v_1},\\ \frac{1}{V_1\left(\mathrm{\infty }\right)}\frac{\mathrm{\partial }{V}_1\left(\mathrm{\infty }\right)}{\mathrm{\partial }{v}_1}& =-{\sigma }_1\frac{1}{S_1\left(\mathrm{\infty }\right)}\frac{N_1{S}_1\left(\mathrm{\infty }\right)}{S_1\left(0\right)}+\frac{{\sigma }_1}{1-v_1}+\frac{1}{v_1},\\ \frac{\mathrm{\partial }{V}_1\left(\mathrm{\infty }\right)}{\mathrm{\partial }{v}_1}& =\frac{V_1\left(\mathrm{\infty }\right)}{v_1}>0.\end{array}$ This is the case that ${\sigma }_1=1$, which vaccine has the lowest efficacy. If the vaccine efficacy is upgraded, people get better protection, so $\frac{\mathrm{\partial }{V}_1\left(\mathrm{\infty }\right)}{\mathrm{\partial }{v}_1}>0,$ holds for all $0<\sigma <1$.

If ${\sigma }_1=0$, which is the case that the vaccine is perfect effective, then (33) $\begin{array}{rl}& \frac{N_1}{S_2\left(\mathrm{\infty }\right)}\left[(c_2{c}_3-c_1{c}_4)\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)\left(\frac{{\sigma }_1{V}_1\left(\mathrm{\infty }\right)}{S_1\left(0\right)}+\frac{V_1\left(\mathrm{\infty }\right)}{V_1\left(0\right)}\right)\right.\\ & \left.+c_1\left(\frac{{\sigma }_1{V}_1\left(\mathrm{\infty }\right)}{S_1\left(0\right)}+\frac{V_1\left(\mathrm{\infty }\right)}{V_1\left(0\right)}\right)+\frac{c_4}{S_1\left(0\right)}\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)-\frac{1}{S_1\left(0\right)}\right]\\ & =\frac{N_1}{S_2\left(\mathrm{\infty }\right)S_1\left(0\right)}\left[(c_2{c}_3-c_1{c}_4)\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)\left({\sigma }_1{V}_1\left(\mathrm{\infty }\right)+\frac{V_1\left(\mathrm{\infty }\right)S_1\left(0\right)}{V_1\left(0\right)}\right)\right.\\ & \left.+c_1({\sigma }_1{V}_1\left(\mathrm{\infty }\right)+\frac{V_1\left(\mathrm{\infty }\right)S_1\left(0\right)}{V_1\left(0\right)})+c_4\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)-1\right]\\ & =\frac{N_1}{S_2\left(\mathrm{\infty }\right)S_1\left(0\right)}\left[(c_2{c}_3-c_1{c}_4)\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)({\sigma }_1{V}_1\left(0\right)+S_1\left(0\right))\right.\\ & \left.+c_1({\sigma }_1{V}_1\left(0\right)+S_1\left(0\right))+c_4\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)-1\right].\end{array}$(33) If ${\sigma }_2$ is also close to 0, then Equation (33(33) $\begin{array}{rl}& \frac{N_1}{S_2\left(\mathrm{\infty }\right)}\left[(c_2{c}_3-c_1{c}_4)\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)\left(\frac{{\sigma }_1{V}_1\left(\mathrm{\infty }\right)}{S_1\left(0\right)}+\frac{V_1\left(\mathrm{\infty }\right)}{V_1\left(0\right)}\right)\right.\\ & \left.+c_1\left(\frac{{\sigma }_1{V}_1\left(\mathrm{\infty }\right)}{S_1\left(0\right)}+\frac{V_1\left(\mathrm{\infty }\right)}{V_1\left(0\right)}\right)+\frac{c_4}{S_1\left(0\right)}\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)-\frac{1}{S_1\left(0\right)}\right]\\ & =\frac{N_1}{S_2\left(\mathrm{\infty }\right)S_1\left(0\right)}\left[(c_2{c}_3-c_1{c}_4)\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)\left({\sigma }_1{V}_1\left(\mathrm{\infty }\right)+\frac{V_1\left(\mathrm{\infty }\right)S_1\left(0\right)}{V_1\left(0\right)}\right)\right.\\ & \left.+c_1({\sigma }_1{V}_1\left(\mathrm{\infty }\right)+\frac{V_1\left(\mathrm{\infty }\right)S_1\left(0\right)}{V_1\left(0\right)})+c_4\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)-1\right]\\ & =\frac{N_1}{S_2\left(\mathrm{\infty }\right)S_1\left(0\right)}\left[(c_2{c}_3-c_1{c}_4)\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)({\sigma }_1{V}_1\left(0\right)+S_1\left(0\right))\right.\\ & \left.+c_1({\sigma }_1{V}_1\left(0\right)+S_1\left(0\right))+c_4\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)-1\right].\end{array}$(33) ) $\begin{array}{rl}& \approx \frac{S_1}{S_2\left(\mathrm{\infty }\right)N_1\left(0\right)}\left[\right(c_2{c}_3-c_1{c}_4\left)\right(S_2\left(0\right)+{\sigma }_2{V}_2\left(0\right)\left)({\sigma }_1{V}_1\right(0)+S_1(0\left)\right)\\ & +c_1({\sigma }_1{V}_1\left(0\right)+S_1\left(0\right))+c_4\left(S_2\right(0)+{\sigma }_2{V}_2(0\left)\right)-1].\end{array}$ From the previous lemma, we know that $\frac{\mathrm{\partial }{S}_1\left(\mathrm{\infty }\right)}{\mathrm{\partial }{v}_1}>0.$

We establish the conclusion.

Theorem 4.1.

In the heterogeneous SVIR model with imperfect vaccine, if we keep the vaccine coverage of subgroup 2 in a steady level, and only change the vaccine coverage level of subgroup 1. Then, we have that when the vaccine has low effectivity to subgroup 1, if vaccine coverage level of subgroup 1 increases, $S_1\left(\mathrm{\infty }\right)$ will decrease, and $V_1\left(\mathrm{\infty }\right)$ will increase; both $S_2\left(\mathrm{\infty }\right)$ and $V_2\left(\mathrm{\infty }\right)$ increase. When the vaccine has high effectivity to subgroup 1, if vaccine coverage level of subgroup 1 increases, $S_1\left(\mathrm{\infty }\right),$ $V_1\left(\mathrm{\infty }\right),$ $S_2\left(\mathrm{\infty }\right)$ and $V_2\left(\mathrm{\infty }\right)$ will increase.

Remark 4.1.

We notice that in Equation (30(30) $\begin{array}{rl}& \frac{N_1}{S_2\left(\mathrm{\infty }\right)}\left[(c_2{c}_3-c_1{c}_4)\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)\left(\frac{{\sigma }_1{V}_1\left(\mathrm{\infty }\right)}{S_1\left(0\right)}+\frac{V_1\left(\mathrm{\infty }\right)}{V_1\left(0\right)}\right)\right.\\ & \left.+c_1\left(\frac{{\sigma }_1{V}_1\left(\mathrm{\infty }\right)}{S_1\left(0\right)}+\frac{V_1\left(\mathrm{\infty }\right)}{V_1\left(0\right)}\right)+\frac{c_4}{S_1\left(0\right)}\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)-\frac{1}{S_1\left(0\right)}\right].\end{array}$(30) ), if (34) ${\sigma }_1{V}_1\left(\mathrm{\infty }\right)+\frac{V_1\left(\mathrm{\infty }\right)S_1\left(0\right)}{V_1\left(0\right)}=\frac{1-c_4\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)}{(c_2{c}_3-c_1{c}_4)\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)+c_1},$(34) there is the change of dynamic behaviour of $S_1\left(\mathrm{\infty }\right)$, we call Equation (34(34) ${\sigma }_1{V}_1\left(\mathrm{\infty }\right)+\frac{V_1\left(\mathrm{\infty }\right)S_1\left(0\right)}{V_1\left(0\right)}=\frac{1-c_4\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)}{(c_2{c}_3-c_1{c}_4)\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)+c_1},$(34) ) the critical-relation equation. We know that both ${\sigma }_1{V}_1\left(\mathrm{\infty }\right)+V_1\left(\mathrm{\infty }\right)S_1\left(0\right)/V_1\left(0\right)$ and $1-c_4\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)$ are positive. If $c_2{c}_3-c_1{c}_4<0$, $(c_2{c}_3-c_1{c}_4)\left(S_2\right(\mathrm{\infty })+{\sigma }_2{V}_2(\mathrm{\infty }\left)\right)+c_1$ may be negative. In such a case, the critical-relation equation does not exit, which means that the relation between ${\sigma }_1$ and ${\sigma }_2$ does not affect the change of $S_1\left(\mathrm{\infty }\right)$, only the value of ${\sigma }_1$ matters. This conclusion is reasonable, because if $c_2{c}_3-c_1{c}_4<0$ holds, we have $c_1$ and $c_4$ much bigger than $c_2$ and $c_3$, that means two subgroups have few contacts with each other. To the extent, we have two separate homogeneous mixing population. If $c_2{c}_3-c_1{c}_4$ is close to 0, or $c_2{c}_3-c_1{c}_4\ge 0$, the critical-relation equation may exist, that means two different subgroups have frequent contact, then the relation of ${\sigma }_1$ and ${\sigma }_2$ will affect the results.

For attack ratios, deriving from the final size relation, we have the conclusion.

Theorem 4.2.

In the heterogeneous SVIR model, if the vaccine coverage level in subgroup 1 increases, all attack ratios decrease; for the increase of vaccine coverage level in subgroup 2, the results hold. The infection rates ${\sigma }_1$ and ${\sigma }_2$ are unrelated to the change of attack ratios.

Remark 4.2.

Herd immunity threshold is not considered in the analysis. In heterogeneous mixing population, the herd immunity threshold will be a relation which the vaccine coverage levels must satisfy, instead of a value.

4.3. Case of perfect vaccine

For the special case of perfect vaccine, we hold all other assumptions. By employing the same technique, we obtain the final size relations (35) $\begin{array}{rl}\ln \frac{S_1\left(0\right)}{S_1\left(\mathrm{\infty }\right)}& =\frac{a_1{p}_{11}}{N_1}\frac{(1-v_1)N_1-S_1\left(\mathrm{\infty }\right)}{{\alpha }_1}+\frac{a_1{p}_{12}}{N_2}\frac{(1-v_2)N_2-S_2\left(\mathrm{\infty }\right)}{{\alpha }_2},\\ \ln \frac{S_2\left(0\right)}{S_2\left(\mathrm{\infty }\right)}& =\frac{a_2{p}_{21}}{N_1}\frac{(1-v_1)N_1-S_1\left(\mathrm{\infty }\right)}{{\alpha }_1}+\frac{a_2{p}_{22}}{N_2}\frac{(1-v_2)N_2-S_2\left(\mathrm{\infty }\right)}{{\alpha }_2},\end{array}$(35) with $v_1$ and $v_2$ representing the vaccine level coverage in subgroup 1 and subgroup 2, respectively. By analysing the final size relations, we have similar results.

Corollary 4.1.

In the heterogeneous SIR model, if the vaccine coverage level in any subgroup increases, $S_{\mathrm{\infty }}$ in both subgroups increase.

Corollary 4.2.

In the heterogeneous SIR model, if the vaccine coverage level in any subgroup increases, the attack ratios in every subgroup also decrease.

5. Extension in general compartmental models

For general compartmental models which may include Exposed class, Quarantine class, Treatment class, etc., we use the general age-of-infection models. By proving the similarity of final size relations obtained from the age-of infection models and from the SIR or SVIR models in homogeneous mixing population and heterogeneous mixing population, we can extend our results in previous sections to the general compartmental models. The same as in [9], we consider a population with a distribution of subpopulation identified by the variable ζ ranging over a ‘state’ space Ω, having sizes $N(t,\zeta )$, respectively, at time t. Suppose that each group ζ member makes $a\left(\zeta \right)$ contacts sufficient to transmit infection per unit time. We assume that the fraction of contacts made by a member of group ζ with a member of group η is $p(\zeta ,\eta )$, then: ${\int }_{\mathrm{\Omega }}p(\zeta ,\eta )\mathrm{d}\eta =1,$ for each ζ. The total number of contacts made in unit time by member of group ζ with members of group η is $a\left(\zeta \right)p(\zeta ,\eta )N(t,\zeta ),$ and by balance relation we have $a\left(\zeta \right)p(\zeta ,\eta )N(t,\zeta )=a\left(\eta \right)p(\eta ,\zeta )N(t,\eta ),$ with $\zeta \ne \eta$.

We assume that there is no permanent movement between groups and that there are no disease deaths, so that $N(t,\zeta )$ is a constant function $N\left(\zeta \right)$ of t for all ζ. We consider two cases, which are perfect vaccine and imperfect vaccine.

1. A multi-group age of infection model with general mixing would be (36) $\begin{array}{rl}{S}^{\mathrm{\prime }}(t,\zeta )& =-a\left(\zeta \right)S(t,\zeta ){\int }_{\mathrm{\Omega }}p(\zeta ,\eta )\frac{\phi (t,\eta )}{N\left(\eta \right)}\mathrm{d}\eta ,\\ \phi (t,\zeta )& ={\phi }_0(t,\zeta )+{\int }_0^{\mathrm{\infty }}[-S^{\mathrm{\prime }}(t-\tau ,\zeta \left)\right]A(\tau ,\zeta )\mathrm{d}\tau .\end{array}$(36) The initial conditions are $S_0\left(\zeta \right)=(1-p(\zeta \left)\right)N\left(\zeta \right)$ and the vaccine is perfectly effective.

   It is possible to obtain a system of final size relations, which is (37) $\begin{array}{rl}& \ln \frac{S(0,\zeta )}{S(\mathrm{\infty },\zeta )}\\ & =a\left(\zeta \right){\int }_{\mathrm{\Omega }}\frac{p(\zeta ,\eta )}{N\left(\eta \right)}\left({\int }_0^{\mathrm{\infty }}{\phi }_0(t,\eta )\mathrm{d}t+\left(S\right(0,\eta )-S(\mathrm{\infty },\eta \left)\right){\int }_0^{\mathrm{\infty }}A\left(\eta \right)\mathrm{d}t\right)\mathrm{d}\eta .\end{array}$(37)

2. We now investigate the case which the vaccine are imperfect, and we model this by including factors $\sigma \left(\zeta \right)$, where $0<\sigma \left(\zeta \right)<1$. The multi-group age of infection model with general mixing (38) $\begin{array}{rl}{S}^{\mathrm{\prime }}(t,\zeta )& =-a\left(\zeta \right)S(t,\zeta ){\int }_{\mathrm{\Omega }}p(\zeta ,\eta )\frac{\phi (t,\eta )}{N\left(\eta \right)}\mathrm{d}\eta ,\\ V^{\mathrm{\prime }}(t,\zeta )& =-\sigma \left(\zeta \right)a\left(\zeta \right)V(t,\zeta ){\int }_{\mathrm{\Omega }}p(\zeta ,\eta )\frac{\phi (t,\eta )}{N\left(\eta \right)}\mathrm{d}\eta ,\\ \phi (t,\zeta )& ={\phi }_0(t,\zeta )+{\int }_0^{\mathrm{\infty }}[-S^{\mathrm{\prime }}(t-\tau ,\zeta -V^{\mathrm{\prime }}(t-\tau ,\zeta )\left]A\right(\tau ,\zeta )\mathrm{d}\tau .\end{array}$(38)

   The final size relations of the model (38(38) $\begin{array}{rl}{S}^{\mathrm{\prime }}(t,\zeta )& =-a\left(\zeta \right)S(t,\zeta ){\int }_{\mathrm{\Omega }}p(\zeta ,\eta )\frac{\phi (t,\eta )}{N\left(\eta \right)}\mathrm{d}\eta ,\\ V^{\mathrm{\prime }}(t,\zeta )& =-\sigma \left(\zeta \right)a\left(\zeta \right)V(t,\zeta ){\int }_{\mathrm{\Omega }}p(\zeta ,\eta )\frac{\phi (t,\eta )}{N\left(\eta \right)}\mathrm{d}\eta ,\\ \phi (t,\zeta )& ={\phi }_0(t,\zeta )+{\int }_0^{\mathrm{\infty }}[-S^{\mathrm{\prime }}(t-\tau ,\zeta -V^{\mathrm{\prime }}(t-\tau ,\zeta )\left]A\right(\tau ,\zeta )\mathrm{d}\tau .\end{array}$(38) ) are (39) $\begin{array}{rl}& \ln \frac{S(0,\zeta )}{S(\mathrm{\infty },\zeta )}=a\left(\zeta \right){\int }_{\mathrm{\Omega }}\frac{p(\zeta ,\eta )}{N\left(\eta \right)}\\ & \times \left({\int }_0^{\mathrm{\infty }}{\phi }_0(t,\eta )\mathrm{d}t+\left(N\right(\eta )-S(\mathrm{\infty },\eta )-V(\mathrm{\infty },\eta \left)\right){\int }_0^{\mathrm{\infty }}A\left(\eta \right)\mathrm{d}t\right)\mathrm{d}\eta ,\\ & \ln \frac{V(0,\zeta )}{V(\mathrm{\infty },\zeta )}=\sigma \left(\zeta \right)a\left(\zeta \right){\int }_{\mathrm{\Omega }}\frac{p(\zeta ,\eta )}{N\left(\eta \right)}\\ & \times \left({\int }_0^{\mathrm{\infty }}{\phi }_0(t,\eta )\mathrm{d}t+\left(N\right(\eta )-S(\mathrm{\infty },\eta )-V(\mathrm{\infty },\eta \left)\right){\int }_0^{\mathrm{\infty }}A\left(\eta \right)\mathrm{d}t\right)\mathrm{d}\eta .\end{array}$(39)

If we compare the final size relations (37(37) $\begin{array}{rl}& \ln \frac{S(0,\zeta )}{S(\mathrm{\infty },\zeta )}\\ & =a\left(\zeta \right){\int }_{\mathrm{\Omega }}\frac{p(\zeta ,\eta )}{N\left(\eta \right)}\left({\int }_0^{\mathrm{\infty }}{\phi }_0(t,\eta )\mathrm{d}t+\left(S\right(0,\eta )-S(\mathrm{\infty },\eta \left)\right){\int }_0^{\mathrm{\infty }}A\left(\eta \right)\mathrm{d}t\right)\mathrm{d}\eta .\end{array}$(37) ) and (39(39) $\begin{array}{rl}& \ln \frac{S(0,\zeta )}{S(\mathrm{\infty },\zeta )}=a\left(\zeta \right){\int }_{\mathrm{\Omega }}\frac{p(\zeta ,\eta )}{N\left(\eta \right)}\\ & \times \left({\int }_0^{\mathrm{\infty }}{\phi }_0(t,\eta )\mathrm{d}t+\left(N\right(\eta )-S(\mathrm{\infty },\eta )-V(\mathrm{\infty },\eta \left)\right){\int }_0^{\mathrm{\infty }}A\left(\eta \right)\mathrm{d}t\right)\mathrm{d}\eta ,\\ & \ln \frac{V(0,\zeta )}{V(\mathrm{\infty },\zeta )}=\sigma \left(\zeta \right)a\left(\zeta \right){\int }_{\mathrm{\Omega }}\frac{p(\zeta ,\eta )}{N\left(\eta \right)}\\ & \times \left({\int }_0^{\mathrm{\infty }}{\phi }_0(t,\eta )\mathrm{d}t+\left(N\right(\eta )-S(\mathrm{\infty },\eta )-V(\mathrm{\infty },\eta \left)\right){\int }_0^{\mathrm{\infty }}A\left(\eta \right)\mathrm{d}t\right)\mathrm{d}\eta .\end{array}$(39) ) with (35(35) $\begin{array}{rl}\ln \frac{S_1\left(0\right)}{S_1\left(\mathrm{\infty }\right)}& =\frac{a_1{p}_{11}}{N_1}\frac{(1-v_1)N_1-S_1\left(\mathrm{\infty }\right)}{{\alpha }_1}+\frac{a_1{p}_{12}}{N_2}\frac{(1-v_2)N_2-S_2\left(\mathrm{\infty }\right)}{{\alpha }_2},\\ \ln \frac{S_2\left(0\right)}{S_2\left(\mathrm{\infty }\right)}& =\frac{a_2{p}_{21}}{N_1}\frac{(1-v_1)N_1-S_1\left(\mathrm{\infty }\right)}{{\alpha }_1}+\frac{a_2{p}_{22}}{N_2}\frac{(1-v_2)N_2-S_2\left(\mathrm{\infty }\right)}{{\alpha }_2},\end{array}$(35) ) and (20(20) $\begin{array}{rl}\ln \frac{S_1\left(0\right)}{S_1\left(\mathrm{\infty }\right)}& =\frac{a_1{p}_{11}}{{\alpha }_1{N}_1}[N_1-S_1(\mathrm{\infty })-V_1(\mathrm{\infty }\left)\right]+\frac{a_1{p}_{12}}{{\alpha }_2{N}_2}[N_2-S_2(\mathrm{\infty })-V_2(\mathrm{\infty }\left)\right],\\ \ln \frac{V_1\left(0\right)}{V_1\left(\mathrm{\infty }\right)}& ={\sigma }_1\left[\frac{a_1{p}_{11}}{{\alpha }_1{N}_1}[N_1-S_1(\mathrm{\infty })-V_1(\mathrm{\infty }\left)\right]+\frac{a_1{p}_{12}}{{\alpha }_2{N}_2}[N_2-S_2(\mathrm{\infty })-V_2(\mathrm{\infty }\left)\right]\right],\\ \ln \frac{S_2\left(0\right)}{S_2\left(\mathrm{\infty }\right)}& =\frac{a_2{p}_{21}}{{\alpha }_1{N}_1}[N_1-S_1(\mathrm{\infty })-V_1(\mathrm{\infty }\left)\right]+\frac{a_2{p}_{22}}{{\alpha }_2{N}_2}[N_2-S_2(\mathrm{\infty })-V_2(\mathrm{\infty }\left)\right],\\ \ln \frac{V_2\left(0\right)}{V_2\left(\mathrm{\infty }\right)}& ={\sigma }_2\left[\frac{a_2{p}_{21}}{{\alpha }_1{N}_1}[N_1-S_1(\mathrm{\infty })-V_1(\mathrm{\infty }\left)\right]+\frac{a_2{p}_{22}}{{\alpha }_2{N}_2}[N_2-S_2(\mathrm{\infty })-V_2(\mathrm{\infty }\left)\right]\right].\end{array}$(20) ), it is clear that the forms of final size relations are similar. The final size relation to the homogeneous mixing age-of-infection model is straightforward to obtain. The generality of the final size relation formulas guarantees that we can easily extend our results from previous sections to the general compartmental models. More details of similar work can be found in [17]. As an example, in [5], it has been proved that the standard susceptible-exposed-infected-recovered (SEIR) epidemic model can be written in the form of the age-of-infection model. So, the property of attack ratios can be applied to the SEIR model. Furthermore, it can be argued that the vaccination game of epidemic with exposed stage has a unique Nash equilibrium. By a similar technique, vaccination games of more complicated epidemics have unique Nash equilibria as well.

6. Conclusion and future work

We demonstrated two vaccination games in homogeneous mixing population and find the Nash equilibrium in each game. To prove that the Nash equilibrium in each game is unique, it is important to consider the monotonicity of the attack ratio. We have formulated the deterministic compartmental vaccination models, in homogeneous mixing population and in heterogeneous population. The final size relations can be obtained from the models. By analysing the final size relations, we proved that the attack ratios can be expressed in terms of the vaccine coverage levels, and also proved that the attack ratios are decreasing functions of the vaccine coverage levels. In other words, we have provided the mathematical proof that the vaccine can protect both susceptible population and vaccinated population better, this is the so-called direct and indirect protection stemming from vaccine. The probability of every individual to survive through the epidemic is greater by increasing the vaccine coverage level. We propose some possible extensions:

• (A1) Compare the measures of derivatives to conduct the sensitivity analysis. If one derivative is bigger than the other, from the standpoint of public health management, to vaccinate one certain subgroup is the priority, since the same amount of vaccination can lower the attack ratio more.

• (A2) The analysis in Section 4 only derives analytical results for the vaccine efficacy factor which is close to 1 and close to 0. Because of too many parameters in the model, there are no explicit formulas to be given. By other methods, we may obtain more accurate results, one possible way is to consider proportional mixing or preferred mixing, instead of general mixing.

• (A3) The vaccination game in heterogeneous mixing population is much more complicated to analyse. The uniqueness of Nash equilibrium does not only depend on the monotonicity of attack ratios, but also depend on the mixing patterns of the population. This part of work needs more effort.

• (A4) We prove that in heterogeneous mixing population, if two subgroups are mainly separate, the vaccine efficacy for subgroup 1 does not affect subgroup 2, vice versa.

Acknowledgments

This work was first inspired by Professor Fred Brauer. The author thanks Fred for weekly discussions on this subject. The author also thanks one anonymous reviewer for helpful comments. This research was supported by NSERC grant.

Disclosure statement

No potential conflict of interest was reported by the authors.