Impact of Imperfect Vaccine, Vaccine Trade-Off and Population Turnover on Infectious Disease Dynamics

Abstract

Vaccination is an essential tool for the management of infectious diseases. However, many vaccines are imperfect, having only a partial protective effect in decreasing disease transmission and/or favouring recovery of infected individuals and possibly exhibiting a trade-off between these two properties. Furthermore, the success of vaccination also depends on the population turnover, and the rate of entry to and exit from the population. We here investigate by means of a mathematical model the interplay between these factors to predict optimal vaccination strategies. We first compute the basic reproduction number and study the global stability of the equilibria. We then assess the most influential parameters determining the total number of infected over time using a sensitivity analysis. We derive conditions for the vaccination coverage and efficiency to achieve disease eradication, assuming different intensities of population turnover (weak and strong), vaccine properties (transmission and/or recovery) and the trade-off between the latter. We show that the minimum vaccination coverage increases with lower population turnover decreases with higher vaccine efficiency (transmission or recovery) and is increased/decreased by up to 15% depending on the vaccine trade-off. We conclude that the coverage target for vaccination campaigns should be evaluated based on the interplay between these factors.

1. Introduction

Vaccination is one of the most effective public health policies for protecting humans and animals from infectious diseases. Global vaccination campaigns have helped eradicate diseases such as smallpox, measles, poliomyelitis and rinderpest in most parts of the world, ultimately saving the lives of millions of humans and animals. By definition, a perfect vaccine would keep vaccinated individuals from becoming infected when exposed to the pathogen. An imperfect vaccine, however, does not prevent vaccinated individuals from becoming infected upon pathogen exposure but may still be beneficial in various ways [1]. For example, imperfect vaccines may provide benefits such as preventing infection, limiting parasite within-host growth and thus reducing the damage done to the host [2] or preventing transmission by infected hosts [3]. As we have seen recently with the epidemic of COVID-19, imperfect vaccines can be used to reduce the number of infected individuals and also to protect individuals at risk of developing a more lethal form of the infection. The use of imperfect vaccines may be advantageous when the vaccination efficiency is volatile and decreases due to the appearance of new variants of the virus [4,5,6].

The effectiveness of a given vaccine is determined not only by its biochemical and immunological properties but also by how the vaccine is deployed and what other health management (biosecurity) measures are in place. Maintaining herd immunity during a disease outbreak, for example, has been promoted as a highly effective disease control strategy [7,8,9]. However, a continuous influx of new susceptible, possibly unvaccinated individuals contributes to the long-term persistence of the disease in the population [10,11]. The frequent introduction of pathogens into a partially immune population (with an intermediate level of population immunity) can lead to longer-lasting epidemics and/or a higher total number of infectious individuals than the introduction into a naive population [10]. This phenomenon is named “epidemic enhancement” [10]. More generally, the population turnover rate, that is, the rate at which individuals can enter and exit the considered population, may affect the effectiveness of control strategies [12,13,14]. In humans and also domesticated animals, population turnover takes the form of immigration and emigration in and out of the population, as well as the birth and death of individuals. The turnover is an often neglected factor in epidemiology when generalising predictions of disease modelling from human to domesticated and wild animal populations.

Moreover, a second parameter of importance in studying the efficiency of vaccination strategies is the existence of biological trade-offs in the epidemiology of infectious diseases. The prime example is the trade-off between parasite virulence and transmission rate, which raises challenges for vaccine manufacturing. Indeed, in the seminal paper by [3], vaccines affecting disease transmission are predicted to possibly lead to a decrease in parasite virulence, while other types of vaccines (reducing within-host growth rate) may lead to an increase in parasite virulence and thus the counter-effect of a worst epidemiological outcome. Interestingly, much work has been devoted to generating precise predictions for virulence evolution in known parasite species by incorporating empirical characterisations of vaccine effects into models capturing the epidemiological details of a given system [15,16,17]. In contrast, the biochemical and immunological trade-offs of the vaccine itself have received little attention. We specifically mean here that vaccination can affect several aspects of the disease dynamics, such as within-host growth and transmission, with possible trade-offs between these characteristics. For example, a vaccine reducing within-host growth may be more or less effective in reducing disease transmission. We, therefore, consider the definition of imperfect vaccines as (i) providing partial protection (non-maximal efficiency) against infection (decreasing transmission) and (ii) partially enhancing the recovery of infected individuals. We are interested in the possible trade-off between these two properties. There has been remarkably little work performed to generally assess how the interplay between different vaccine properties, trade-offs and vaccination strategies influences the burden of the epidemic in a heterogeneous community/population with imperfect vaccination.

The aim of this study is, therefore, to assess, through mathematical modelling, whether the use of vaccines decreasing the infection rate is more efficient to eradicate the disease in a heterogeneous community than a vaccine that both reduces the infection and favours recovery or a vaccine reducing the infection rate but favouring recovery. We also want to assess whether these results depend on the effect of population turnover in order to generalise our results to animal populations. The paper is organised as follows. First, the model is formulated in Section 2. We then compute the basic properties of the steady-state solutions as well as the existence of local and global stability of the equilibria of the model (Section 3). We then perform a numerical sensitivity analysis of the model and study examples of numerical analyses for different parameter values to describe the interaction between population turnover and vaccine trade-offs on the epidemiological outcome. We conclude by providing predictions on the applicability of these results to vaccination strategies in human populations, but also domesticated (and wild) animal species for which turnover rates represent a different end of a continuum.

2. Model Formulation

The formulation of the model is based on compartmental modelling [18], which consists of creating virtual reservoirs called compartments. A compartment is a kinetically homogeneous structure. This means that any individual who enters a compartment is identical, from the epidemiological point of view to any other already present in that compartment. A mathematical model, therefore, consists of describing the flow of individuals between the various compartments.

To study the dynamics of an infectious disease during and after the vaccination campaign, we modify the model formulated in [3] by adding a recovered compartment, and we consider a frequency-dependent disease transmission (incidence rate). The model takes into account only host-to-host transmission of the disease. Since many vaccines do not guarantee perfect immunity, we consider a heterogeneous host community/population with two types of hosts: fully susceptible to the disease or partially resistant to infection due to the imperfect vaccination. The fully susceptible hosts consist of uninfected ($S_1$) and infected ($I_1$) individuals. Among the partially resistant hosts, there are uninfected ($S_2$) and infected ($I_2$) individuals. All infected individuals (fully susceptible or partially resistant) can become recovered (R), and all recovered individuals are fully immune to reinfection [19]. Thus, the total population at time t, $N\left(t\right)$ is given by

$$N\left(t\right)=S_1\left(t\right)+S_2\left(t\right)+I_1\left(t\right)+I_2\left(t\right)+R\left(t\right).$$

We assume the parasite population to be monomorphic (having only one type or genotype). We also assume that new uninfected hosts arise through birth and immigration at a constant rate, $\theta \ge 0$. Among these new uninfected, a proportion, $0\le p\le 1$, is partially immune due to the vaccination, while the remaining proportion, $1-p$, is susceptible (completely vulnerable to the parasite). Uninfected, infected and recovered hosts die naturally at a rate $\mu \ge 0$, and infected hosts suffer additional mortality due to the virulence of the parasite. Since host resistance due to vaccination may reduce the impact of the parasite within-host growth [3], we assume the virulence of the parasite on fully susceptible hosts, $d_1\ge 0$, to be greater than that on partially resistant hosts, $d_2\ge 0$. Uninfected hosts become infected with the forces of infection ${\lambda }_1\left(t\right)={\beta }_{11}{\displaystyle \frac{I_1\left(t\right)}{N\left(t\right)}}+{\beta }_{12}{\displaystyle \frac{I_2\left(t\right)}{N\left(t\right)}}$ and ${\lambda }_2\left(t\right)={\beta }_{21}{\displaystyle \frac{I_1\left(t\right)}{N\left(t\right)}}+{\beta }_{22}{\displaystyle \frac{I_2\left(t\right)}{N\left(t\right)}}$ when they are fully susceptible or partially resistant, respectively. The rates of transmission are ${\beta }_{11}\ge 0$ (${\beta }_{21}\ge 0$) from infected, $I_1$, to susceptible individuals $S_1$ ($S_2$), while ${\beta }_{12}\ge 0$ (${\beta }_{22}\ge 0$) is the transmission rate from infected, $I_2$, to susceptible individuals $S_1$ ($S_2$). Since the resistance can decrease the probability of becoming infected [3], we generally assume ${\beta }_{21}\le {\beta }_{11}$ and ${\beta }_{22}\le {\beta }_{12}$. Recovery rates may differ between the fully susceptible ${\gamma }_1\ge 0$ and the partially resistant host, ${\gamma }_2\ge 0$. The schematic diagram of the model is shown in Figure 1.

Mathematically, the model is as follows:

$$\left\{\begin{array}{ccc}\hfill {\displaystyle \frac{\mathrm{d}{S}_1}{\mathrm{d}t}}\hfill & \hfill =\hfill & \theta (1-p)-{\lambda }_1\left(t\right)S_1\left(t\right)-\mu S_1\left(t\right),\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}{S}_2}{\mathrm{d}t}}\hfill & \hfill =\hfill & \theta p-{\lambda }_2\left(t\right)S_2\left(t\right)-\mu S_2\left(t\right),\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}\hfill & \hfill =\hfill & {\lambda }_1\left(t\right)S_1\left(t\right)-(\mu +{\gamma }_1+d_1)I_1\left(t\right),\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}\hfill & \hfill =\hfill & {\lambda }_2\left(t\right)S_2\left(t\right)-(\mu +{\gamma }_2+d_2)I_2\left(t\right),\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}R}{\mathrm{d}t}}\hfill & \hfill =\hfill & {\gamma }_1{I}_1\left(t\right)+{\gamma }_2{I}_2\left(t\right)-\mu R\left(t\right).\hfill \end{array}\right.$$

(1)

A summary of the biological significance of the model parameters (1) is given in

Table 1. Description and value of the model parameters.

Parameter	Description	Units	Value	Source
$\theta$	Recruitment rate	person·day${}^{-1}$	variable	Assumed
$\mu$	Natural mortality rate	day${}^{-1}$	variable	Assumed
p	Proportion of new hosts vaccinated	-	variable	Assumed
${\beta }_{11}$	Transmission rate from $I_1$ to $S_1$	day${}^{-1}$	variable	Assumed
${\beta }_{12}$	Transmission rate from $I_2$ to $S_1$	day${}^{-1}$	variable	Assumed
${\beta }_{21}$	Transmission rate from $I_1$ to $S_2$	day${}^{-1}$	variable	Assumed
${\beta }_{22}$	Transmission rate from $I_2$ to $S_2$	day${}^{-1}$	variable	Assumed
$d_1$	Mortality rate due to infection of $S_1$	day${}^{-1}$	$0.0008$	[9]
$d_2$	Mortality rate due to infection of $S_2$	day${}^{-1}$	$0.0001$	[9]
${\gamma }_1$	Recovery rate of $I_1$	day${}^{-1}$	$0.1$	[9]
${\gamma }_2$	Recovery rate of $I_2$	day${}^{-1}$	$0.13$	[9]

.

3. Mathematical Analysis

3.1. Basic Properties

First, we study the basic characteristics of the system solutions: the existence, non-negativity and boundedness of solutions. These are (1) essential to make sure that the model (1) is well defined mathematically and epidemiologically and (2) useful for the proofs of the stability results.

3.1.1. Positive Invariance of the Non-Negative Orthant

For any associated Cauchy problem, system (1), which is a $C^{\infty }$-differentiable system, has a unique maximal solution.

Lemma 1. 

The following result corresponds to Proposition B.7, Appendix B in [20]. Let D be an open subset of ${\mathbb{R}}^n$, $f:\mathbb{R}\times D\to {\mathbb{R}}^n$, be a vector-valued function, $f=(f_1,f_2,\cdots ,f_n).$ Consider a system of ODEs of the form

$$x^{\prime }=f(t,x).$$

(2)

Suppose that f in Equation (2) has the property that solutions of initial value problems $x\left(t_0\right)=x_0\ge 0$ are unique and for all i$f_i(t,x)\ge 0$ whenever $x\ge 0$ satisfies $x_i=0$. Then $x\left(t\right)\ge 0$ for all $t\ge t_0$ for which it is defined, provided $x\left(t_0\right)\ge 0.$

Theorem 1. 

If the initial conditions of system (1) are such that $S_1\left(0\right)\ge 0$, $S_2\left(0\right)\ge 0$, $I_1\left(0\right)\ge 0$, $I_2\left(0\right)\ge 0$ and $R\left(0\right)\ge 0$, then the solution $(S_1\left(t\right),S_2\left(t\right),I_1\left(t\right),I_2\left(t\right),R\left(t\right))$ of the system equation is non-negative for all $t\ge 0.$

Proof. 

Considering model (1). We have

$$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{d}{S}_1}{\mathrm{d}t}}{|}_{S_1=0}=& \theta (1-p)\ge 0,\hfill \\ \\ \hfill {\displaystyle \frac{\mathrm{d}{S}_2}{\mathrm{d}t}}{|}_{S_2=0}=& \theta p\ge 0,\hfill \\ \\ \hfill {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}{|}_{I_1=0}=& {\beta }_{11}{\displaystyle \frac{I_1\left(t\right)}{N\left(t\right)}}{S}_1\left(t\right)\ge 0,\hfill \\ \\ \hfill {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}{|}_{I_2=0}=& {\beta }_{21}{\displaystyle \frac{I_1\left(t\right)}{N\left(t\right)}}{S}_2\left(t\right)\ge 0,\hfill \\ \\ \hfill {\displaystyle \frac{\mathrm{d}R}{\mathrm{d}t}}{|}_{R=0}=& {\gamma }_1{I}_1\left(t\right)+{\gamma }_2{I}_2\left(t\right)\ge 0,\hfill \end{array}$$

for all $S_1,S_2,I_1,I_2,R\ge 0$. By using Lemma 1, we conclude that the solution $(S_1\left(t\right),S_2\left(t\right),I_1\left(t\right),I_2\left(t\right),R\left(t\right))$ of the system equation is non-negative for all $t\ge 0.$ □

Thus, solutions of system (1) with non-negative initial conditions will be non-negative for all $t\ge 0$.

3.1.2. Boundedness of Solutions

Since the variables of model (1) are non-negative, and we are dealing with the dynamic of a number of individuals, it is important and biologically realistic that the total number of individuals does not explode (that is, it is bounded).

Lemma 2. 

The closed set

$$\mathsf{\Omega }=\left\{(S_1\left(t\right),S_2\left(t\right),I_1\left(t\right),I_2\left(t\right),R\left(t\right))\in {\mathbb{R}}^5:S_1\left(t\right)\ge 0,S_2\left(t\right)\ge 0,I_1\left(t\right)\ge 0,I_2\left(t\right)\ge 0,R\left(t\right)\ge 0,N\left(t\right)⩽{\displaystyle \frac{\theta }{\mu }}\right\}$$

is positively invariant and attracting for system (1).

Proof. 

Using system (1), the dynamics of the total human population satisfies:

$${\displaystyle \frac{\mathrm{d}N}{\mathrm{d}t}}=\theta -\mu N-d_1{I}_1-d_2{I}_2⩽\theta -\mu N.$$

Integrating both sides of the expression above, we deduce that

$$N\left(t\right)⩽{\displaystyle \frac{\theta }{\mu }}+\left(N\left(0\right)-{\displaystyle \frac{\theta }{\mu }}\right)e^{-\mu t},\forall t⩾0,$$

(3)

where $N\left(0\right)$ is the value of $N\left(t\right)$ at time zero. We deduce that if $N\left(0\right)⩽{\displaystyle \frac{\theta }{\mu }}$, then $0⩽N\left(t\right)⩽{\displaystyle \frac{\theta }{\mu }}$, $\forall t⩾0$ and $\mathsf{\Omega }$ is positively invariant. If $N\left(0\right)⩾{\displaystyle \frac{\theta }{\mu }}$, then from Equation (3) the total population decreases, and the solutions of system (1) enter $\mathsf{\Omega }$. Hence $N\left(t\right)$ is bounded as $t\to \infty$, which means that $\mathsf{\Omega }$ is attracting. □

Remark 1. 

We know from Theorem 13 in [21] that every maximal solution of the Cauchy problem (2) that is bounded is global; that is, it exists for all $t\ge 0$. Then, every maximal solution of system (1) is well-defined for all $t\ge 0$.

System (1) is epidemiologically and mathematically well-posed in $\mathsf{\Omega }$ since its state variables are non-negative, and the size of the total population is bounded. The maximum value of N represents the size of the total population under the ideal situation without infection.

3.2. Disease-Free Equilibrium and Its Stability

For the analysis of the spread of an infection, we define the disease-free equilibrium (DFE), which is the state of the population without infection. The disease-free equilibrium is deduced from the resolution of system (1) by taking $I_1=0$ and $I_2=0$. Thus, the disease-free equilibrium satisfies the following system of equations:

$$\left\{\begin{array}{c}\theta (1-p)-\mu S_1^0=0,\hfill \\ \theta p-\mu S_2^0=0.\hfill \end{array}\right.$$

(4)

Solving the system of equations (4) yields the disease-free equilibrium point:

$$Q^0=(S_1^0,S_2^0,0,0,0),$$

where $S_1^0={\displaystyle \frac{\theta (1-p)}{\mu }}$, $S_2^0={\displaystyle \frac{\theta p}{\mu }}$ and $N^0=S_1^0+S_2^0={\displaystyle \frac{\theta }{\mu }}$.

The linear stability of $Q^0$ depends on the well-known reproduction number ${\mathcal{R}}_0$, which is defined as the average number of secondary cases caused by an infected individual during its infectious period when introduced into a population of susceptible individuals. We study the stability of the equilibrium through the next-generation operator [22,23]. Recalling the notations in [23] for system (1), the matrices $\mathcal{F}$ of the new infection and $\mathcal{V}$ of the remaining transfer terms are, respectively, given by

$$\mathcal{F}=\left[\begin{array}{c}{\beta }_{11}{\displaystyle \frac{S_1{I}_1}{N}}+{\beta }_{12}{\displaystyle \frac{S_1{I}_2}{N}}\\ \\ {\beta }_{21}{\displaystyle \frac{S_2{I}_1}{N}}+{\beta }_{22}{\displaystyle \frac{S_2{I}_2}{N}}\end{array}\right]\mathrm{and}\mathcal{V}=\left[\begin{array}{c}(\mu +{\gamma }_1+d_1)I_1\\ \\ (\mu +{\gamma }_2+d_2)I_2\end{array}\right].$$

The Jacobian matrices of $\mathcal{F}$ and $\mathcal{V}$ at $Q^0$ are, respectively,

$$F=\left[\begin{array}{cc}{\beta }_{11}{\displaystyle \frac{S_1^0}{N^0}}& {\beta }_{12}{\displaystyle \frac{S_1^0}{N^0}}\\ \\ {\beta }_{21}{\displaystyle \frac{S_2^0}{N^0}}& {\beta }_{22}{\displaystyle \frac{S_2^0}{N^0}}\end{array}\right]\mathrm{and}V=\left[\begin{array}{cc}\mu +{\gamma }_1+d_1& 0\\ \\ 0& \mu +{\gamma }_2+d_2\end{array}\right].$$

(5)

Then,

$$F{V}^{-1}=\left[\begin{array}{cc}{\displaystyle \frac{{\beta }_{11}{S}_1^0}{N^0(\mu +{\gamma }_1+d_1)}}& {\displaystyle \frac{{\beta }_{12}{S}_1^0}{N^0(\mu +{\gamma }_2+d_2)}}\\ \\ {\displaystyle \frac{{\beta }_{21}{S}_2^0}{N^0(\mu +{\gamma }_1+d_1)}}& {\displaystyle \frac{{\beta }_{22}{S}_2^0}{N^0(\mu +{\gamma }_2+d_2)}}\end{array}\right],$$

and the reproduction number of model system (1) is

$$\begin{array}{cc}\hfill {\mathcal{R}}_0=\rho \left(F{V}^{-1}\right)=& \hfill {\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{S_1^0}{N^0}}{\mathcal{R}}_{0,11}+{\displaystyle \frac{S_2^0}{N^0}}{\mathcal{R}}_{0,22}+\sqrt{{\left({\displaystyle \frac{S_1^0}{N^0}}{\mathcal{R}}_{0,11}-{\displaystyle \frac{S_2^0}{N^0}}{\mathcal{R}}_{0,22}\right)}^2+4{\displaystyle \frac{S_1^0}{N^0}}{\displaystyle \frac{S_2^0}{N^0}}{\mathcal{R}}_{0,12}{\mathcal{R}}_{0,21}}\right],\hfill \\ \hfill {\mathcal{R}}_0=& \hfill {\displaystyle \frac{1}{2}}\left[(1-p){\mathcal{R}}_{0,11}+p{\mathcal{R}}_{0,22}+\sqrt{{\left((1-p){\mathcal{R}}_{0,11}-p{\mathcal{R}}_{0,22}\right)}^2+4p(1-p){\mathcal{R}}_{0,12}{\mathcal{R}}_{0,21}}\right],\hfill \end{array}$$

(6)

where ${\displaystyle \frac{S_1^0}{N^0}}=1-p$ (${\displaystyle \frac{S_2^0}{N^0}}=p$) is the proportion of susceptible individuals that have not been vaccinated (have been vaccinated) at the DFE $Q^0$. Similarly, we define ${\mathcal{R}}_{0,11}={\displaystyle \frac{{\beta }_{11}}{\mu +{\gamma }_1+d_1}}$ as the average number of secondary cases generated by an unvaccinated infected individual during its infectious period through interaction with the unvaccinated population. Furthermore, ${\mathcal{R}}_{0,12}={\displaystyle \frac{{\beta }_{12}}{\mu +{\gamma }_1+d_1}}$ represents the average number of secondary cases generated by a vaccinated infected individual in the unvaccinated part of the population, ${\mathcal{R}}_{0,21}={\displaystyle \frac{{\beta }_{21}}{\mu +{\gamma }_2+d_2}}$ is the average number of secondary cases generated by an unvaccinated infected individual in the vaccinated part of the population, and ${\mathcal{R}}_{0,22}={\displaystyle \frac{{\beta }_{22}}{\mu +{\gamma }_2+d_2}}$ represents the average number of secondary cases generated by an infected vaccinated individual in the vaccinated part of the population. Further, $\rho \left(F{V}^{-1}\right)$ is the spectral radius of $F{V}^{-1}$.

Remark 2. 

From the expression of the reproduction number ${\mathcal{R}}_0$ in Equation (6), we deduce that ${\mathcal{R}}_0\ge max\{(1-p){\mathcal{R}}_{0,11};p{\mathcal{R}}_{0,22}\}$. Moreover, using (6) for $p=0$ (all new hosts are unvaccinated), ${\mathcal{R}}_0={\mathcal{R}}_{0,11}$. Further, if $p=1$ (all new hosts are vaccinated), then ${\mathcal{R}}_0={\mathcal{R}}_{0,22}$.

The importance of the reproduction number is due to the result given in the next lemma derived from Theorem 2 in [23].

Lemma 3. 

The disease-free equilibrium $Q^0$ of system (1) is locally asymptotically stable in Ω if ${\mathcal{R}}_0<1$, and unstable if ${\mathcal{R}}_0>1$.

The biological meaning of Lemma 3 is that a sufficiently small number of infected hosts does not induce an epidemic unless the reproduction number ${\mathcal{R}}_0$ is greater than unity. That is, the disease rapidly dies out (when ${\mathcal{R}}_0<1$) if the initial number of infected hosts is in the basin of attraction of the DFE, $Q^0$. Global asymptotic stability of the DFE is required to better control the disease. In addition, analysing the expansion of the basin of attraction of $Q^0$ is a more challenging task for the model under consideration, involving a fairly new result. For this purpose, we use Theorems 2.1 and 2.2 in [24].

Theorem 2. 

If ${\mathcal{R}}_0⩽1$, the disease-free equilibrium $Q^0$ of system (1) is globally asymptotic stable in Ω. If ${\mathcal{R}}_0>1$, $Q^0$ is unstable, system (1) is uniformly persistent, and there exists at least one endemic equilibrium in the interior of Ω.

Proof. 

See Appendix A. □

As a consequence of the meaning of Theorem 2 and Remark 2, we can confidently deduce that the disease can be eradicated from the host community if the value of ${\mathcal{R}}_0$ is reduced to less than unity, independently of whether individuals introduced to the population are all vaccinated or not.

3.3. Endemic Equilibrium and Its Stability

Let $Q^{*}=(S_1^{*},S_2^{*},I_1^{*},I_2^{*},R^{*})$ be the positive endemic equilibrium (EE) of model system (1). Then, the positive endemic equilibrium can be obtained by setting the right-hand side of all equations in model system (1) to zero, giving:

$$\left\{\begin{array}{c}\theta (1-p)-{\beta }_{11}{\displaystyle \frac{S_1^{*}{I}_1^{*}}{N^{*}}}-{\beta }_{12}{\displaystyle \frac{S_1^{*}{I}_2^{*}}{N^{*}}}-\mu S_1^{*}=0,\hfill \\ \\ \theta p-{\beta }_{21}{\displaystyle \frac{S_2^{*}{I}_1^{*}}{N^{*}}}-{\beta }_{22}{\displaystyle \frac{S_2^{*}{I}_2^{*}}{N^{*}}}-\mu S_2^{*}=0,\hfill \\ \\ {\beta }_{11}{\displaystyle \frac{S_1^{*}{I}_1^{*}}{N^{*}}}+{\beta }_{12}{\displaystyle \frac{S_1^{*}{I}_2^{*}}{N^{*}}}-(\mu +{\gamma }_1+d_1)I_1^{*}=0,\hfill \\ \\ {\beta }_{21}{\displaystyle \frac{S_2^{*}{I}_1^{*}}{N^{*}}}+{\beta }_{22}{\displaystyle \frac{S_2^{*}{I}_2^{*}}{N^{*}}}-(\mu +{\gamma }_2+d_2)I_2^{*}=0,\hfill \\ \\ {\gamma }_1{I}_1^{*}+{\gamma }_2{I}_2^{*}-\mu R^{*}=0.\hfill \end{array}\right.$$

(7)

Given the complexity of system (7), we are not determining an explicit formula for the endemic equilibrium point $Q^{*}$. Note that determining $Q^{*}$ is often very difficult to be carried out when the system is complex and has a large size. However, to prove the existence of $Q^{*}$, we can rewrite system (7) as a fixed point problem and use Theorem 2.1 in [25]. To do this, we solve system (7). After algebraic manipulations, we obtain: $R^{*}={\displaystyle \frac{{\gamma }_1{I}_1^{*}+{\gamma }_2{I}_2^{*}}{\mu }}$, $S_1^{*}={\displaystyle \frac{\theta (1-p)N^{*}}{{\beta }_{11}{I}_1^{*}+{\beta }_{12}{I}_2^{*}+\mu N^{*}}}$, $S_2^{*}={\displaystyle \frac{\theta p{N}^{*}}{{\beta }_{21}{I}_1^{*}+{\beta }_{22}{I}_2^{*}-d_1{I}_1^{*}-d_2{I}_2^{*}+\theta }}$, $I_1^{*}={\displaystyle \frac{\theta (1-p)({\beta }_{11}{I}_1^{*}+{\beta }_{12}{I}_2^{*})}{(\mu +{\gamma }_1+d_1)({\beta }_{11}{I}_1^{*}+{\beta }_{12}{I}_2^{*}-d_1{I}_1^{*}-d_2{I}_2^{*}+\theta )}}=H_1\left(I^{*}\right)$ and $I_2^{*}={\displaystyle \frac{\theta p({\beta }_{21}{I}_1^{*}+{\beta }_{22}{I}_2^{*})}{(\mu +{\gamma }_2+d_2)({\beta }_{21}{I}_1^{*}+{\beta }_{22}{I}_2^{*}-d_1{I}_1^{*}-d_2{I}_2^{*}+\theta )}}=H_2\left(I^{*}\right)$ with $I^{*}=(I_1^{*},I_2^{*}).$

Then, the endemic equilibrium is the fixed points of H given by $I=H\left(I\right)$ where $I=(I_1,I_2)$. By definition, H is continuous, monotonously non-decreasing and strictly sublinear. H is also a bounded function that maps the non-negative orthant $\mathsf{\Omega }$ into itself. Moreover, $H\left(0\right)=0$ by definition and the jacobian of H at the zero, $H^{{}^{\prime }}\left(0\right)$, exists and is irreducible since

$$H^{{}^{\prime }}\left(0\right)=\left[\begin{array}{cc}{\beta }_{11}{a}_1& {\beta }_{12}{a}_1\\ \\ {\beta }_{21}{a}_2& {\beta }_{22}{a}_2\end{array}\right]=F{V}^{-1},$$

where $a_1={\displaystyle \frac{1-p}{\mu +{\gamma }_1+d_1}}$ and $a_2={\displaystyle \frac{p}{\mu +{\gamma }_2+d_2}}.$

We deduce that the spectral radius $\rho \left(H^{{}^{\prime }}\left(0\right)\right)$ of the matrix $H^{{}^{\prime }}\left(0\right)$ is ${\mathcal{R}}_0.$ Then, the existence and the uniqueness of a non-negative fixed point occur if and only if ${\mathcal{R}}_0>1$.

Proposition 1. 

System (1) has only one endemic equilibrium whenever ${\mathcal{R}}_0>1$.

We establish the following result to analyse the stability of $Q^{*}$.

Theorem 3. 

If ${\mathcal{R}}_0>1$, the endemic equilibrium $Q^{*}$ is globally asymptotic stable in Ω.

Proof. 

Consider the following Lyapunov candidate function:

$$L=L_1+L_2+L_3+L_4,$$

where $L_1=S_1-S_1^{*}-S_1^{*}log\left({\displaystyle \frac{S_1}{S_1^{*}}}\right)$, $L_2=S_2-S_2^{*}-S_2^{*}log\left({\displaystyle \frac{S_2}{S_2^{*}}}\right)$, $L_3=I_1-I_1^{*}-I_1^{*}log\left({\displaystyle \frac{I_1}{I_1^{*}}}\right)$ and $L_4=I_2-I_2^{*}-I_2^{*}log\left({\displaystyle \frac{I_2}{I_2^{*}}}\right)$.

Using the inequality $1-z+log\left(z\right)⩽0$ for $z>0$ with equality if and only if $z=1$, differentiation and using the EE values give

$$L^{{}^{\prime }}=L_1^{{}^{\prime }}+L_2^{{}^{\prime }}+L_3^{{}^{\prime }}+L_4^{{}^{\prime }},$$

where

$$\begin{array}{cc}\hfill L_1^{{}^{\prime }}\hfill & =\left(1-{\displaystyle \frac{S_1^{*}}{S_1}}\right){\displaystyle \frac{\mathrm{d}{S}_1}{\mathrm{d}t}}\hfill \\ \\ \hfill & =\left(1-{\displaystyle \frac{S_1^{*}}{S_1}}\right)\left[{\beta }_{11}{\displaystyle \frac{S_1^{*}{I}_1^{*}}{N^{*}}}-{\beta }_{11}{\displaystyle \frac{S_1{I}_1}{N}}+{\beta }_{12}{\displaystyle \frac{S_1^{*}{I}_2^{*}}{N^{*}}}-{\beta }_{12}{\displaystyle \frac{S_1{I}_2}{N}}-\mu S_1+\mu S_1^{*}\right]\hfill \\ \\ \hfill & =-{\displaystyle \frac{\mu {(S_1-S_1^{*})}^2}{S_1}}+{\beta }_{11}{\displaystyle \frac{S_1^{*}{I}_1^{*}}{N^{*}}}\left[1-{\displaystyle \frac{S_1^{*}}{S_1}}-{\displaystyle \frac{S_1{I}_1{N}^{*}}{S_1^{*}{I}_1^{*}N}}+{\displaystyle \frac{I_1{N}^{*}}{I_1^{*}N}}\right]+{\beta }_{12}{\displaystyle \frac{S_1^{*}{I}_2^{*}}{N^{*}}}\left[1-{\displaystyle \frac{S_1^{*}}{S_1}}-{\displaystyle \frac{S_1{I}_2{N}^{*}}{S_1^{*}{I}_2^{*}N}}+{\displaystyle \frac{I_2{N}^{*}}{I_2^{*}N}}\right].\hfill \\ \\ \hfill \mathrm{Then}{L}_1^{{}^{\prime }}⩽\hfill & {\beta }_{11}{\displaystyle \frac{S_1^{*}{I}_1^{*}}{N^{*}}}\left[{\displaystyle \frac{I_1{N}^{*}}{I_1^{*}N}}-log\left({\displaystyle \frac{I_1{N}^{*}}{I_1^{*}N}}\right)-{\displaystyle \frac{S_1{I}_1{N}^{*}}{S_1^{*}{I}_1^{*}N}}+log\left({\displaystyle \frac{S_1{I}_1{N}^{*}}{S_1^{*}{I}_1^{*}N}}\right)\right]\hfill \\ \\ \hfill \hfill & +{\beta }_{12}{\displaystyle \frac{S_1^{*}{I}_2^{*}}{N^{*}}}\left[{\displaystyle \frac{I_2{N}^{*}}{I_2^{*}N}}-log\left({\displaystyle \frac{I_2{N}^{*}}{I_2^{*}N}}\right)-{\displaystyle \frac{S_1{I}_2{N}^{*}}{S_1^{*}{I}_2^{*}N}}+log\left({\displaystyle \frac{S_1{I}_2{N}^{*}}{S_1^{*}{I}_2^{*}N}}\right)\right].\hfill \end{array}$$

(8)

We can also deduce in an analogous way:

$$\begin{array}{cc}\hfill L_2^{{}^{\prime }}\hfill & \hfill ⩽{\beta }_{22}{\displaystyle \frac{S_2^{*}{I}_2^{*}}{N^{*}}}\left[{\displaystyle \frac{I_2{N}^{*}}{I_2^{*}N}}-log\left({\displaystyle \frac{I_2{N}^{*}}{I_2^{*}N}}\right)-{\displaystyle \frac{S_2{I}_2{N}^{*}}{S_2^{*}{I}_2^{*}N}}+log\left({\displaystyle \frac{S_2{I}_2{N}^{*}}{S_2^{*}{I}_2^{*}N}}\right)\right]\hfill \\ \\ \hfill \hfill & \hfill +{\beta }_{21}{\displaystyle \frac{S_2^{*}{I}_1^{*}}{N^{*}}}\left[{\displaystyle \frac{I_1{N}^{*}}{I_1^{*}N}}-log\left({\displaystyle \frac{I_1{N}^{*}}{I_1^{*}N}}\right)-{\displaystyle \frac{S_2{I}_1{N}^{*}}{S_2^{*}{I}_1^{*}N}}+log\left({\displaystyle \frac{S_2{I}_1{N}^{*}}{S_2^{*}{I}_1^{*}N}}\right)\right].\hfill \end{array}$$

(9)

We also have

$$\begin{array}{cc}\hfill L_3^{{}^{\prime }}& =\left(1-{\displaystyle \frac{I_1^{*}}{I_1}}\right){\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}\hfill \\ \\ \hfill \hfill & =\left(1-{\displaystyle \frac{I_1^{*}}{I_1}}\right)\left[{\beta }_{11}{\displaystyle \frac{S_1{I}_1}{N}}+{\beta }_{12}{\displaystyle \frac{S_1{I}_2}{N}}-(\mu +{\gamma }_1+d_1)I_1\right]\hfill \\ \\ \hfill \hfill & =\left(1-{\displaystyle \frac{I_1^{*}}{I_1}}\right)\left[{\beta }_{11}{\displaystyle \frac{S_1{I}_1}{N}}+{\beta }_{12}{\displaystyle \frac{S_1{I}_2}{N}}-{\beta }_{11}{\displaystyle \frac{S_1^{*}{I}_1}{N^{*}}}+{\beta }_{12}{\displaystyle \frac{S_1^{*}{I}_2^{*}{I}_1}{N^{*}{I}_1^{*}}}\right]\hfill \\ \\ \hfill \hfill & ={\beta }_{11}{\displaystyle \frac{S_1^{*}{I}_1^{*}}{N^{*}}}\left[{\displaystyle \frac{S_1{I}_1{N}^{*}}{S_1^{*}{I}_1^{*}N}}-{\displaystyle \frac{S_1{N}^{*}}{S_1^{*}N}}-{\displaystyle \frac{I_1}{I_1^{*}}}+1\right]+{\beta }_{12}{\displaystyle \frac{S_1^{*}{I}_2^{*}}{N^{*}}}\left[{\displaystyle \frac{S_1{I}_2{N}^{*}}{S_1^{*}{I}_2^{*}N}}-{\displaystyle \frac{S_1{I}_1^{*}{I}_2{N}^{*}}{S_1^{*}{I}_1{I}_2^{*}N}}-{\displaystyle \frac{I_1}{I_1^{*}}}+1\right],\hfill \\ \\ \hfill L_3^{{}^{\prime }}& ⩽{\beta }_{11}{\displaystyle \frac{S_1^{*}{I}_1^{*}}{N^{*}}}\left[{\displaystyle \frac{S_1{I}_1{N}^{*}}{S_1^{*}{I}_1^{*}N}}-log\left({\displaystyle \frac{S_1{I}_1{N}^{*}}{S_1^{*}{I}_1^{*}N}}\right)-{\displaystyle \frac{I_1}{I_1^{*}}}+log\left({\displaystyle \frac{I_1}{I_1^{*}}}\right)\right]\hfill \\ \\ \hfill \hfill & +{\beta }_{12}{\displaystyle \frac{S_1^{*}{I}_2^{*}}{N^{*}}}\left[{\displaystyle \frac{S_1{I}_2{N}^{*}}{S_1^{*}{I}_2^{*}N}}-log\left({\displaystyle \frac{S_1{I}_2{N}^{*}}{S_1^{*}{I}_2^{*}N}}\right)-{\displaystyle \frac{I_1}{I_1^{*}}}+log\left({\displaystyle \frac{I_1}{I_1^{*}}}\right)\right].\hfill \end{array}$$

(10)

Similarly, we obtain

$$\begin{array}{cc}\hfill L_4^{{}^{\prime }}\hfill & \hfill ⩽{\beta }_{22}{\displaystyle \frac{S_2^{*}{I}_2^{*}}{N^{*}}}\left[{\displaystyle \frac{S_2{I}_2{N}^{*}}{S_2^{*}{I}_2^{*}N}}-log\left({\displaystyle \frac{S_2{I}_2{N}^{*}}{S_2^{*}{I}_2^{*}N}}\right)-{\displaystyle \frac{I_2}{I_2^{*}}}+log\left({\displaystyle \frac{I_2}{I_2^{*}}}\right)\right]\hfill \\ \\ \hfill \hfill & \hfill +{\beta }_{21}{\displaystyle \frac{S_2^{*}{I}_1^{*}}{N^{*}}}\left[{\displaystyle \frac{S_2{I}_1{N}^{*}}{S_2^{*}{I}_1^{*}N}}-log\left({\displaystyle \frac{S_2{I}_1{N}^{*}}{S_2^{*}{I}_1^{*}N}}\right)-{\displaystyle \frac{I_2}{I_2^{*}}}+log\left({\displaystyle \frac{I_2}{I_2^{*}}}\right)\right].\hfill \end{array}$$

(11)

Therefore, by adding (8)–(11), we deduce

$$\begin{array}{cc}\hfill L^{{}^{\prime }}⩽& \left(-{\displaystyle \frac{I_1{N}^{*}}{I_1^{*}N}}+log\left({\displaystyle \frac{I_1{N}^{*}}{I_1^{*}N}}\right)\right)\left(-{\beta }_{11}{\displaystyle \frac{S_1^{*}{I}_1^{*}}{N^{*}}}-{\beta }_{21}{\displaystyle \frac{S_2^{*}{I}_1^{*}}{N^{*}}}\right)\hfill \\ \hfill & +\left(-{\displaystyle \frac{I_2{N}^{*}}{I_2^{*}N}}+log\left({\displaystyle \frac{I_2{N}^{*}}{I_2^{*}N}}\right)\right)\left(-{\beta }_{12}{\displaystyle \frac{S_1^{*}{I}_2^{*}}{N^{*}}}-{\beta }_{22}{\displaystyle \frac{S_2^{*}{I}_2^{*}}{N^{*}}}\right)\hfill \\ \hfill & +\left(-{\displaystyle \frac{I_1}{I_1^{*}}}+log\left({\displaystyle \frac{I_1}{I_1^{*}}}\right)\right)\left({\beta }_{11}{\displaystyle \frac{S_1^{*}{I}_1^{*}}{N^{*}}}+{\beta }_{12}{\displaystyle \frac{S_1^{*}{I}_2^{*}}{N^{*}}}\right)\hfill \\ \hfill & +\left(-{\displaystyle \frac{I_2}{I_2^{*}}}+log\left({\displaystyle \frac{I_2}{I_2^{*}}}\right)\right)\left({\beta }_{22}{\displaystyle \frac{S_2^{*}{I}_2^{*}}{N^{*}}}+{\beta }_{21}{\displaystyle \frac{S_2^{*}{I}_1^{*}}{N^{*}}}\right).\hfill \\ \hfill \mathrm{Then}{L}^{{}^{\prime }}⩽& 0,\sin \mathrm{ce}-z+log\left(z\right)⩽-1,\forall z>0.\hfill \end{array}$$

Since $\left\{Q^{*}\right\}$ is the only invariant subset in $\mathsf{\Omega }$ where $L=0$, therefore, by La Salle’s invariance principle [26], $Q^{*}$ is globally global asymptotic stable in $\mathsf{\Omega }$. □

The epidemiological consequence of this theorem is that the disease persists as endemic in the host population as soon as ${\mathcal{R}}_0>1$.

3.4. Herd Immunity Threshold

Herd immunity is a form of indirect protection from infectious disease that occurs when a sufficient percentage of a population has become immune to an infection, whether through previous infections or vaccination and thereby reducing the likelihood of infection for individuals lacking immunity. This is due to the fact that immune individuals are unlikely to contribute to disease transmission, disrupting chains of infection, which stops or slows down the spread of disease. To compute the herd immunity threshold associated with system (1), we set the reproduction number, ${\mathcal{R}}_0$, to one and solve for $p={\displaystyle \frac{S_2^0}{N^0}}$, which is the proportion of susceptible individuals that have been vaccinated at the DFE, $Q^0.$ Then we have,

$$\begin{array}{cc}\hfill {\mathcal{R}}_0=1⟺& {\left[2-{\mathcal{R}}_{0,11}+({\mathcal{R}}_{0,11}-{\mathcal{R}}_{0,22})p\right]}^2={\left[{\mathcal{R}}_{0,11}-({\mathcal{R}}_{0,11}+{\mathcal{R}}_{0,11})p\right]}^2+4p(1-p){\mathcal{R}}_{0,12}{\mathcal{R}}_{0,21}\hfill \\ \hfill ⟺& \left[{({\mathcal{R}}_{0,11}-{\mathcal{R}}_{0,22})}^2-{({\mathcal{R}}_{0,11}+{\mathcal{R}}_{0,22})}^2+4{\mathcal{R}}_{0,12}{\mathcal{R}}_{0,21}\right]p^2+[2(2-{\mathcal{R}}_{0,11})({\mathcal{R}}_{0,11}-{\mathcal{R}}_{0,22})\hfill \\ \hfill & +2{\mathcal{R}}_{0,11}({\mathcal{R}}_{0,11}+{\mathcal{R}}_{0,22})-4{\mathcal{R}}_{0,12}{\mathcal{R}}_{0,21}]p+{(2-{\mathcal{R}}_{0,11})}^2-{\mathcal{R}}_{0,11}^2=0.\hfill \end{array}$$

Thus, solving ${\mathcal{R}}_0=1$ is equivalent to finding the roots of polynomial $Q\left(p\right)$ given by:

$$Q\left(p\right)=A{p}^2+Bp+C,$$

(12)

where $A=4{\mathcal{R}}_{0,12}{\mathcal{R}}_{0,21}-4{\mathcal{R}}_{0,11}{\mathcal{R}}_{0,22}$, $B=4{\mathcal{R}}_{0,11}(1+{\mathcal{R}}_{0,22})-4({\mathcal{R}}_{0,22}+{\mathcal{R}}_{0,12}{\mathcal{R}}_{0,21})$ and

$C=4(1-{\mathcal{R}}_{0,11})$.

Noting that negative thresholds are biologically meaningless (in our case), the conditions for $Q\left(p\right)$ to have positive real roots are determined below. For this purpose, we perform a case analysis to determine the positive real zeros of Q.

Let $\Delta =B^2-4AC$ be the discriminant of the equation $Q\left(p\right)=0$.

Case 1 

Suppose $A=0$. Then

$$p_c=-\frac{C}{B}$$

is the only real root of Q. In addition, $p_c>0$ if and only if B and C have opposite signs and $B\ne 0.$

Case 2 

Suppose $A\ne 0$ and $\Delta =0$. Then

$$p_{c_0}=-{\displaystyle \frac{B}{2A}}$$

is the only real root of Q. Further $p_{c_0}>0$ if and only if A and B have opposite signs.

Case 3 

Suppose $A\ne 0$ and $\Delta >0$. Then

$$p_{c_1}={\displaystyle \frac{-B-\sqrt{\Delta }}{2A}}\mathrm{and}{p}_{c_2}={\displaystyle \frac{-B+\sqrt{\Delta }}{2A}}$$

are the real roots of Q.

Moreover, if $A>0$, then

$$\left\{\begin{array}{c}{p}_{c_1}>0\mathrm{if}\mathrm{and}\mathrm{only}\mathrm{if}\sqrt{\Delta }<-B,\hfill \\ p_{c_2}>0\mathrm{if}\mathrm{and}\mathrm{only}\mathrm{if}\sqrt{\Delta }>B.\hfill \end{array}\right.$$

Therefore, Q has two positive real roots if $A>0$, $B<0$, $C>0$ and $\Delta >0$. In addition, it has one positive real root if ($A>0$, $B<0$, $C<0$ and $\Delta >0$) or ($A>0$, $B>0$ and $C<0$ and $\Delta >0$).

Finally, if $A<0$, then

$$\left\{\begin{array}{c}{p}_{c_1}>0\mathrm{if}\mathrm{and}\mathrm{only}\mathrm{if}\sqrt{\Delta }>-B,\hfill \\ p_{c_2}>0\mathrm{if}\mathrm{and}\mathrm{only}\mathrm{if}\sqrt{\Delta }<B.\hfill \end{array}\right.$$

Therefore, Q has two positive real roots if $A<0$, $B>0$, $C<0$ and $\Delta >0$. It has one positive real root if ($A<0$, $B>0$, $C>0$ and $\Delta >0$) or ($A<0$, $B<0$, $C>0$ and $\Delta >0$).

Theorems 2 and 3 can be combined to give the following result:

Corollary 1. 

An imperfect vaccine can lead to the elimination of the disease if $Q\left(p\right)>0$ (i.e., ${\mathcal{R}}_0<1$). If $Q\left(p\right)<0$ (i.e., ${\mathcal{R}}_0>1$), then the disease persists in the population.

The implication of Corollary 1 is that the use of an imperfect vaccine can lead to the elimination of the disease in the host population if the proportion of vaccinated individuals satisfies one of the following conditions:

1.

$p>p_c$, if $A=0$, $B>0$ and $C<0$;

2.

$p\in [0,p_c[$, if $A=0$, $B>0$ and $C>0$;

3.

$p\ne p_{c_0}$, if $A>0$, $\Delta =0$ and $B<0$;

4.

$p\in [0,p_{c_1}[$ or $p>p_{c_2}$, if $A>0$, $\Delta >0$, $B<0$ and $C>0$;

5.

$p>p_{c_1}$ or $p>p_{c_2}$, if ($A>0$, $\Delta >0$, $B<0$ and $C<0$) or ($A>0$, $\Delta >0$, $B>0$ and $C<0$);

6.

$p\in ]p_{c_2},p_{c_1}[$, if $A<0$, $\Delta >0$, $B>0$ and $C<0$;

7.

$p\in [0,p_{c_1}[$ or $p\in [0,p_{c_2}[$, if ($A<0$, $\Delta >0$, $B>0$ and $C>0$) or ($A<0$, $\Delta >0$, $B<0$ and $C>0$).

Conversely, the disease persists in the population if the proportion of individuals vaccinated satisfies one of these conditions:

1.

$p\in [0,p_c[$, if $A=0$, $B>0$ and $C<0$;

2.

$p>p_c$, if $A=0$, $B>0$ and $C>0$;

3.

$p\ne p_{c_0}$, if $A<0$, $\Delta =0$ and $B>0$;

4.

$p\in ]p_{c_1},p_{c_2}[$, if $A>0$, $\Delta >0$, $B<0$ and $C>0$;

5.

$p\in [0,p_{c_1}[$ or $p\in [0,p_{c_2}[$, if ($A>0$, $\Delta >0$, $B<0$ and $C<0$) or ($A>0$, $\Delta >0$, $B>0$ and $C<0$);

6.

$p\in [0,p_{c_2}[$ or $p>p_{c_1}$, if $A<0$, $\Delta >0$, $B>0$ and $C<0$;

7.

$p>p_{c_1}$ or $p>p_{c_2}$, if ($A<0$, $\Delta >0$, $B>0$ and $C>0$) or ($A<0$, $\Delta >0$, $B<0$ and $C>0$).

We conclude the analytical part of our study by stating that the eradication of a disease is conditioned by the proportion of vaccinated individuals; this threshold for vaccination coverage is called the critical vaccination proportion ($p_c$). In some cases, there is one critical proportion that determines whether the basic reproduction number,${\mathcal{R}}_0$, is less than one or not. In other cases, two critical proportions are found defining the occurrence of three different possible dynamics: disease eradication when ${\mathcal{R}}_0<1$, endemic disease dynamics when ${\mathcal{R}}_0>1$ with the presence or absence of epidemiological oscillations of the number of infected individuals. In the latter case of the two thresholds, the analytical results derived above do not allow the prediction of the epidemiological dynamics and the vaccination proportions. We, therefore, provide numerical simulations in the follow-up section.

4. Numerical Simulations

We refine the above analytical results by numerical simulations to assess the influence of the various model parameters and the impact of population turnover and trade-offs in vaccination efficiency on the epidemiological dynamics (i.e., the number of infected individuals and ${\mathcal{R}}_0$). To illustrate the behaviour of model (1), we use parameter values for the mortality rates, $d_1$ and $d_2$, and the recovery rates, ${\gamma }_1,{\gamma }_2$, measured for COVID-19 as an example of a highly transmissible disease (based on data from the United States [9]). In order to assess the influence of the various parameters of the model on the epidemiological outcome, we vary their values as described in Table 1. Note that we do not attempt here to precisely model the COVID-19 epidemic, but we focus on highly transmissible diseases relevant to public health. We indeed aim to go beyond applicability to a particular diseases (COVID-19) and to provide a generalised overview of the influence of vaccination trade-offs on epidemics.

4.1. Global Sensitivity Analysis

Uncertainty/sensitivity analyses are first used to determine which model input parameters have the greatest impact on the epidemiological outcome [27]. The sensitivity analysis of the model parameters is carried out to measure the correlation between the model parameters and (1) the total number of infected individuals ($I_1+I_2$) and (2) the threshold parameter ${\mathcal{R}}_0$. The analysis is performed using the Latin Hypercube Sampling (LHS) technique and partial rank correlation coefficients (PRCCs) [27]. In our analysis, 1000 model simulations are performed by running the model for 200 time steps (equivalent to 200 days), and the number of infected is recorded at time points 50, 100 and 200. To perform the sensitivity analysis, each parameter has a parameter range defined by the maximum (the minimum), being $50\%$ greater (less) than its baseline (values in

Table A1. PRCC of model’s parameters at time t (days) with strong PI and weak efficiency of vaccine. The values $\theta =1000$, $\mu =0.09$, ${\beta }_{11}=0.35$, ${\beta }_{12}=0.28$, ${\beta }_{21}=0.175$, ${\beta }_{22}=0.14$ are used as baseline.

Parameters	Range of Parameters			Total Infected: ${\mathit{I}}_1+{\mathit{I}}_2$
	Min	Baseline	Max	$\mathit{t}=50$ Days	$\mathit{t}=100$ Days	$\mathit{t}=200$ Days
$\theta$	500	1000	1500	$0.71395$ **	$0.78511$ **	$0.76166$ ***
p	0	$0.5$	1	$0.020314$	$0.0029584$	$0.028397$
${\beta }_{11}$	$0.175$	$0.35$	$0.525$	$0.85757$ ***	$0.8731$ ***	$0.87175$ ***
${\beta }_{12}$	$0.14$	$0.28$	$0.42$	$0.0047432$	$0.027724$	$-0.030496$
${\beta }_{21}$	$0.0875$	$0.175$	$0.2625$	$0.0090246$	$-0.012341$	$0.026579$
${\beta }_{22}$	$0.07$	$0.14$	$0.21$	$-0.047262$	$0.02905$	$-0.037461$
$\mu$	$0.045$	$0.09$	$0.135$	$-0.7695$ **	$-0.80652$ ***	$-0.79222$ **
$d_1$	$0.0004$	$0.0008$	$0.0012$	$-0.012188$	$0.03368$	$-0.046922$
$d_2$	$0.00005$	$0.0001$	$0.00015$	$-0.025215$	$0.016188$	$-0.043869$
${\gamma }_1$	$0.05$	$0.1$	$0.15$	$-0.78315$ **	$-0.84015$ ***	$-0.82903$ ***
${\gamma }_2$	$0.0625$	$0.13$	$0.1925$	$0.010702$	$0.05007$	$0.012449$

**: PRCC values: 0.7 to 0.79 or −0.7 to −0.79; ***: PRCC values: 0.8 to 0.99 or −0.8 to −0.99.

,

Table A2. PRCC of model’s parameters at time t days with strong PI and strong efficiency of vaccine, when $\theta =1000$, $\mu =0.09$, ${\beta }_{11}=0.35$, ${\beta }_{12}=0.28$, ${\beta }_{21}=0.035$, ${\beta }_{22}=0.028$ as baseline.

Parameters	Range of Parameters			Total Infected: ${\mathit{I}}_1+{\mathit{I}}_2$
	Min	Baseline	Max	$\mathit{t}=50$ Days	$\mathit{t}=100$ Days	$\mathit{t}=200$ Days
$\theta$	500	1000	1500	$0.7381$ **	$0.76387$ **	$0.78486$ **
p	0	$0.5$	1	$0.003905$	$0.0037356$	$0.027257$
${\beta }_{11}$	$0.175$	$0.35$	$0.525$	$0.86469$ ***	$0.87427$ ***	$0.88181$ ***
${\beta }_{12}$	$0.14$	$0.28$	$0.42$	$0.0079816$	$0.033516$	$0.030158$
${\beta }_{21}$	$0.0175$	$0.035$	$0.0525$	$0.0012134$	$0.018087$	$-0.00058438$
${\beta }_{22}$	$0.014$	$0.028$	$0.042$	$0.021841$	$-0.0038364$	$0.018881$
$\mu$	$0.045$	$0.09$	$0.135$	$-0.78849$ **	$-0.80304$ ***	$-0.81535$ ***
$d_1$	$0.0004$	$0.0008$	$0.0012$	$-0.054627$	$0.066816$	$0.019678$
$d_2$	$0.00005$	$0.0001$	$0.00015$	$-0.033227$	$-0.021472$	$-0.028882$
${\gamma }_1$	$0.05$	$0.1$	$0.15$	$-0.80324$ ***	$-0.83346$ ***	$-0.84421$ ***
${\gamma }_2$	$0.0625$	$0.13$	$0.1925$	$-0.0099732$	$-0.02272$	$0.0020891$

**: PRCC values: 0.7 to 0.79 or −0.7 to −0.79; ***: PRCC values: 0.8 to 0.99 or −0.8 to −0.99.

,

Table A3. PRCC of model’s parameters at time t days with weak PI and weak efficiency of vaccine, when $\theta =10$, $\mu =0.0009$, ${\beta }_{11}=0.35$, ${\beta }_{12}=0.28$, ${\beta }_{21}=0.175$, ${\beta }_{22}=0.14$ as baseline.

Parameters	Range of Parameters			Total Infected: ${\mathit{I}}_1+{\mathit{I}}_2$
	Min	Baseline	Max	$\mathit{t}=50$ Days	$\mathit{t}=100$ Days	$\mathit{t}=200$ Days
$\theta$	5	10	15	$0.45737$	$0.37556$	$0.51163$ *
p	0	$0.5$	1	$-0.041574$	$0.028378$	$0.030938$
${\beta }_{11}$	$0.175$	$0.35$	$0.525$	$-0.63334$ *	$-0.23892$	$0.355$
${\beta }_{12}$	$0.14$	$0.28$	$0.42$	$-0.23979$	$-0.24053$	$-0.13989$
${\beta }_{21}$	$0.0875$	$0.175$	$0.2625$	$-0.90072$ ***	$-0.90502$ ***	$-0.80837$ ***
${\beta }_{22}$	$0.07$	$0.14$	$0.21$	$-0.52059$ *	$-0.50519$ *	$-0.30843$
$\mu$	$0.00045$	$0.0009$	$0.00135$	$-0.031697$	$-0.18722$	$-0.15951$
$d_1$	$0.0004$	$0.0008$	$0.0012$	$0.012078$	$-0.038623$	$0.01511$
$d_2$	$0.00005$	$0.0001$	$0.00015$	$0.028409$	$0.0088495$	$0.047733$
${\gamma }_1$	$0.05$	$0.1$	$0.15$	$-0.12428$	$0.81303$ ***	$0.59284$ *
${\gamma }_2$	$0.0625$	$0.13$	$0.1925$	$0.48726$	$0.62754$ *	$0.48082$

*: PRCC values: 0.5 to 0.69 or −0.5 to −0.69; ***: PRCC values: 0.8 to 0.99 or −0.8 to −0.99.

,

Table A4. PRCC of model’s parameters at time t days with weak PI and strong efficiency of vaccine, when $\theta =10$, $\mu =0.0009$, ${\beta }_{11}=0.35$, ${\beta }_{12}=0.28$, ${\beta }_{21}=0.035$, ${\beta }_{22}=0.028$ as baseline.

Parameters	Range of Parameters			Total Infected: ${\mathit{I}}_1+{\mathit{I}}_2$
	Min	Baseline	Max	$\mathit{t}=50$ Days	$\mathit{t}=100$ Days	$\mathit{t}=200$ Days
$\theta$	5	10	15	$0.5751$ *	$0.48818$	$0.61256$ *
p	0	$0.5$	1	$0.052835$	$0.014154$	$0.050557$
${\beta }_{11}$	$0.175$	$0.35$	$0.525$	$-0.70943$ **	$-0.47854$	$0.44458$
${\beta }_{12}$	$0.14$	$0.28$	$0.42$	$-0.16357$	$-0.16371$	$-0.03405$
${\beta }_{21}$	$0.0175$	$0.035$	$0.0525$	$-0.84854$ ***	$-0.90973$ ***	$-0.85731$ ***
${\beta }_{22}$	$0.014$	$0.028$	$0.042$	$-0.17909$	$-0.22613$	$-0.14446$
$\mu$	$0.00045$	$0.0009$	$0.00135$	$-0.40329$	$-0.61646$ *	$-0.72754$ **
$d_1$	$0.0004$	$0.0008$	$0.0012$	$-0.072168$	$0.04039$	$-0.040258$
$d_2$	$0.00005$	$0.0001$	$0.00015$	$0.019586$	$-0.053637$	$-0.030518$
${\gamma }_1$	$0.05$	$0.1$	$0.15$	$-0.81298$ ***	$0.76378$ **	$0.69479$ *
${\gamma }_2$	$0.0625$	$0.13$	$0.1925$	$0.20028$	$0.31528$	$0.2891$

*: PRCC values: 0.5 to 0.69 or −0.5 to −0.69; **: PRCC values: 0.7 to 0.79 or −0.7 to −0.79; ***: PRCC values: 0.8 to 0.99 or −0.8 to −0.99.

). We then divide each parameter range into 1000 equally large sub-intervals, and draw a value per parameter within that interval using a Uniform draw. By this means, we obtain a uniform distribution of 1000 parameter values for each parameter. The parameter space (or LHS matrix) has dimensions of length 11, with each dimension specifying an uncertain parameter vector of length 1000. The base parameter values are chosen to define several scenarios of interest regarding the intensity of the turnover (weak and strong) and the efficiency of the vaccine (weak and strong). In PRCC analysis, the parameters with larger positive or negative PRCC values (>0.5 or <$-0.5$) and with correspondingly small p-values (<0.05) are deemed the most influential in determining the outcome of the model. A positive (negative) correlation coefficient corresponds to an increasing (decreasing) monotonic trend between the chosen response function and the parameter under consideration. The results of the PRCC analyses are found in Table A1, Table A2, Table A3, Table A4 in Appendix B.

Based on the results in Table A1, Table A2, Table A3, Table A4, we provide, in

Table 2. Summary of the influence of parameters on the total numbers of infected at different time points.

Scenarios	Total Infected: ${\mathit{I}}_1+{\mathit{I}}_2$
	$\mathit{t}=50$ Days	$\mathit{t}=100$ Days	$\mathit{t}=200$ Days
Strong turnover and weak efficiency	$\theta$(+), ${\beta }_{11}(+)$, $\mu (-)$, ${\gamma }_1(-)$	$\theta (+)$, ${\beta }_{11}(+)$, $\mu (-)$, ${\gamma }_1(-)$	$\theta (+)$, ${\beta }_{11}(+)$, $\mu (-)$, ${\gamma }_1(-)$
Strong turnover and strong efficiency	$\theta (+)$, ${\beta }_{11}(+)$, $\mu (-)$, ${\gamma }_1(-)$	$\theta (+)$, ${\beta }_{11}(+)$, $\mu (-)$, ${\gamma }_1(-)$	$\theta (+)$, ${\beta }_{11}(+)$, $\mu (-)$, ${\gamma }_1(-)$
Weak turnover and weak efficiency	${\beta }_{11}(-)$, ${\beta }_{21}(-)$, ${\beta }_{22}(-)$	${\beta }_{21}(-)$, ${\beta }_{22}(-)$, ${\gamma }_1(+)$, ${\gamma }_2(+)$	$\theta (+)$, ${\beta }_{21}(-)$, ${\gamma }_1(+)$
Weak turnover and strong efficiency	$\theta (+)$, ${\beta }_{11}(-)$, ${\beta }_{21}(-)$, ${\gamma }_1(-)$	${\beta }_{21}(-)$, $\mu (-)$, ${\gamma }_1(+)$	$\theta (+)$, ${\beta }_{21}(-)$, $\mu (-)$, ${\gamma }_1(+)$

, a summary of the parameters that significantly affect the number of infected individuals. Overall, it appears that the recruitment rate, $\theta$, and the recovery rate of the infected who have not been vaccinated, ${\gamma }_1$, are the two main parameters driving the number of infected individuals. This suggests that an effective control strategy should aim to significantly limit the immigration of new hosts in the population (to decrease $\theta$) and improve the treatment of infected individuals (to increase ${\gamma }_1$). We then proceed to a similar analysis with ${\mathcal{R}}_0$ and summarise the sensitivity analysis of the LHS and PRCC techniques in Figure 2. We find, perhaps unsurprisingly, that the proportion of new hosts vaccinated, p, is the most significant parameter explaining the change in ${\mathcal{R}}_0$, along with the transmission rate from unvaccinated infected to unvaccinated susceptibles, ${\beta }_{11}$, and the recovery rate of the infected who have not been vaccinated, ${\gamma }_1$ (Table 2).

4.2. Interplay between Vaccine Efficiency and Population Turnover

We now study the effect of population turnover and vaccine efficiency on the epidemiological dynamics. Specifically, we use numerical simulations to find the vaccination coverage necessary to eradicate the disease in the community (${\mathcal{R}}_0$ satisfying the Corollary 1) under two population turnover rates (fixing the ratio $\theta /\mu$, we define strong turnover with $\theta =1000$ and $\mu =0.09$, and weak with $\theta =10$ and $\mu =0.0009$), when the efficiency of the vaccine only reduces transmission. The vaccine efficiency is set as weak (${\beta }_{21}=(1-0.5){\beta }_{11}$ and ${\beta }_{22}=(1-0.5){\beta }_{12}$, defining an efficiency of $50\%$) or strong (${\beta }_{21}=(1-0.9){\beta }_{11}$ and ${\beta }_{22}=(1-0.9){\beta }_{12}$, defining an efficiency of $90\%$).

4.2.1. Strong Population Turnover

The epidemiological dynamics in Figure 3b under strong turnover and weak vaccine efficiency (${\mathcal{R}}_0=1.2352$) show that the dynamics reach the endemic disease equilibrium. Furthermore, if p takes a value between 0 and $p_1$ (with $p_1\approx 0.696$), the basic reproduction number is greater than 1, but if p is between $p_1$ and 1, the basic reproduction number is less than 1 (as predicted in the analytical results in Corollary 1). Therefore, to eradicate the disease under strong population turnover and weak efficiency of the vaccine, a minimum vaccination rate is needed and is defined by $p_1$. Under strong turnover and strong efficiency (Figure 3d, with ${\mathcal{R}}_0=0.9808$) the disease becomes extinct. Furthermore, if parameter p is between 0 and $p_2$ with $p_2\approx 0.489$, the basic reproduction number is greater than 1, while for p between $p_2$ and 1, the basic reproduction number is less than 1. Therefore, to eradicate the disease in this context of strong turnover and strong efficiency of the vaccine, there is a need to vaccinate more than $48.9\%$ of the new host individuals.

4.2.2. Weak Population Turnover

To illustrate a weak population turnover, we consider the values $\theta =10$ and $\mu =0.0009$, noting that the ratio of $\theta /\mu$ is the same as for the strong turnover investigated above. Under weak turnover, the epidemiological dynamics exhibit damped oscillations (recurring outbreaks) before stabilising at the endemic state with disease persistence (Figure 4b with ${\mathcal{R}}_0=2.2551$, Figure 4d with ${\mathcal{R}}_0=1.8276$). These oscillations are due to the fact that individuals migrate rapidly in the recovered compartment, and a new outbreak only occurs when a sufficient number of susceptible individuals are available for new recruitment into the population and recovered individuals lose their immunity (so-called waning immunity). This phenomenon was also described in [7,10,11,28], and the effect of turnover and waning immunity is specifically described in [7,10].

With respect to the control of the disease, under weak vaccine efficiency, p can take any value between 0 and 1, and the basic reproduction number is always greater than 1 (Figure 4b with ${\mathcal{R}}_0=2.2551$). In contrast, when vaccine efficiency is strong, three cases occur, Figure 4d (with ${\mathcal{R}}_0=1.8276$). When p has a value between 0 and $p_3$, with $p_3\approx 0.753$, the basic reproduction number is greater than 1, and we observe a damped periodicity of the number of infected individuals converging towards a stable endemic state. When p takes a value between $p_3$ and $p_4$ (with $p_4\approx 0.756$), the basic reproduction number, ${\mathcal{R}}_0$, is greater than 1, but no oscillations are observed. For $p\in [p_4,1]$, the basic reproduction number, ${\mathcal{R}}_0$, is less than 1, and the disease becomes extinct. Note that between $p_3$ and $p_4$, the behaviour can change very finely, but the resolution of our simulations does not allow us to decide on a very precise bound when oscillations occur or not. Therefore, to eradicate the disease in this context of weak population turnover and strong efficiency of the vaccine, high vaccination coverage (more than $75.6\%$ of the new host individuals) is needed. Our results extend those in [29], showing that it is feasible to control the disease by a weakly efficient vaccine acting on disease transmission but that the required vaccination coverage depends on the population turnover. We note that the persistence of an endemic equilibrium is predicted by the condition ${\mathcal{R}}_0>1$, even if damped oscillations in the number of infected individuals occur. In other words, while the population turnover does not factor directly in the analytical expression of ${\mathcal{R}}_0$, it enters only indirectly by affecting the proportion of susceptible individuals available (Equation 6). The simulation results provide examples of the analytical expressions obtained in Equation (12) following Corollary 1.

4.3. Interplay between Types of Vaccines and Population Turnover

We now assume that a vaccine has two potential mechanisms of action on the disease, namely blocking transmission and/or favouring the recovery of infected individuals. We investigate the effect of these vaccine types on the epidemiology depending on the population turnover. Specifically, model (1) is slightly modified to allow for the assessment of the efficiency of the vaccine regarding the probability of being infected and the recovery rate. This is achieved by simply rescaling the parameters as follows:

$$\begin{array}{c}\hfill {\beta }_{21}=(1-\epsilon ){\beta }_{11},{\beta }_{22}=(1-\epsilon ){\beta }_{12},\mathrm{and}{\gamma }_1=(1-\nu ){\gamma }_2,\end{array}$$

(13)

where $0\le \epsilon \le 1$ represents the effect of the vaccine on disease transmission, and $0\le \nu \le 1$ represents the effect of the vaccine on recovery. Substituting the rescaled expressions in Equation (13) into model (1), one deduces that the basic reproduction number in model (1) can be rewritten as:

$$\begin{array}{cc}\hfill {\mathcal{R}}_0=& \hfill {\displaystyle \frac{1}{2}}\left[(1-p){\mathcal{R}}_{0,11}+p{\mathcal{R}}_{0,22}+\sqrt{{\left((1-p){\mathcal{R}}_{0,11}-p{\mathcal{R}}_{0,22}\right)}^2+4p(1-p){\mathcal{R}}_{0,12}{\mathcal{R}}_{0,21}}\right],\hfill \end{array}$$

(14)

with ${\mathcal{R}}_{0,11}={\displaystyle \frac{{\beta }_{11}}{\mu +(1-\nu ){\gamma }_2+d_1}}$, ${\mathcal{R}}_{0,12}={\displaystyle \frac{{\beta }_{12}}{\mu +(1-\nu ){\gamma }_2+d_1}}$, ${\mathcal{R}}_{0,21}={\displaystyle \frac{(1-\epsilon ){\beta }_{11}}{\mu +{\gamma }_2+d_2}}$ and ${\mathcal{R}}_{0,22}={\displaystyle \frac{(1-\epsilon ){\beta }_{12}}{\mu +{\gamma }_2+d_2}}$. Simulations are carried out to assess the interplay of the type of vaccine and the population turnover.

Under a strong population turnover, as expected, the value of the reproduction number decreases as coverage and efficiency of the vaccine on the transmission increase (Figure 5a), and if the vaccine is designed to only decrease the transmission by $80\%$ (i.e., $\epsilon =0.8$), the eradication of the disease in the host population can be achieved (${\mathcal{R}}_0<1$) if at least $70\%$ of the population is vaccinated (Figure 5a). On the other hand, the value of the reproduction number decreases as coverage increases and the efficiency of the vaccine favouring recovery decreases (Figure 5b). With a vaccine designed to enhance recovery by $20\%$ (i.e., $\nu =0.2$), the eradication of the disease in the host population can be achieved (${\mathcal{R}}_0<1$) if at least $68\%$ of the population is vaccinated (Figure 5b). In Figure 5c, we present the effect of the combined efficiency of the vaccine (decreasing transmission and favouring recovery) on the reproduction number at $p=0.5$. The eradication of the disease can be achieved (${\mathcal{R}}_0<1$) if the vaccine has a combined efficiency of at least $85\%$ against infection (and thus transmission) and at least $20\%$ to enhance recovery (for a given vaccination coverage of $p=0.5$). These figures represent subsets of the general results presented in Figure A1, in which ${\mathcal{R}}_0$ is a function of $\epsilon$, $\nu$ and p. The use of a vaccine with a combined efficiency (decreasing transmission and favouring recovery) can be associated with vaccination coverage in order to achieve the elimination of the disease. For example, with vaccination coverage of $20\%$ ($p=0.2$), it is not possible to eliminate the disease no matter the combined efficiency of the vaccine (Figure A2), while at $80\%$ coverage ($p=0.8$), there are several combinations of vaccine types, decreasing transmission and favouring recovery, that can promote disease control (Figure A2).

The above results dramatically change under a weak population turnover. As expected, the value of the reproduction number decreases as coverage and efficiency of the vaccine on the transmission increase (Figure A3a), but higher vaccination coverage is needed compared to the strong population turnover to achieve ${\mathcal{R}}_0<1$. Moreover, it is not possible to eradicate the disease if (1) the vaccine is only efficient in enhancing recovery, no matter the vaccination coverage (Figure A3b), or (2) if the efficiency of the vaccine is combined, but vaccination coverage is $p=0.5$ (Figure A4). The general results of ${\mathcal{R}}_0$ as a function of $\epsilon$, $\nu$ and p demonstrate that under weak population turnover, disease eradication requires a very strong efficiency of the vaccine and high coverage (Figure A5).

4.4. Interplay between Vaccine Efficiency Trade-Off and Population Turnover

Thus far, we have assumed that all parameters of vaccine efficiency can be independently chosen from one another. We study, here, the epidemiological dynamics when there exists a possible (and realistic) trade-off (relationship) between the vaccine efficiency on transmission and on recovery. We assume three possible trade-off curves: convex ($\nu ={\epsilon }^2$), concave ($\nu =\sqrt{\epsilon }$) or linear ($\nu =\epsilon$). Under a strong population turnover, assuming a vaccine of at least $60\%$ efficiency, disease eradication can be achieved (${\mathcal{R}}_0<1$) if the coverage is at least $65\%$ under a convex trade-off (Figure 6a), at least $80\%$ under a concave trade-off (Figure 6b) and at least $75\%$ under a linear trade-off (Figure 6c). Imposing vaccine trade-off, therefore, affects the shape of the ${\mathcal{R}}_0$ curves in Figure 6a–c compared to Figure 5a,b, and may be important to predict the minimum vaccination coverage to be achieved. However, under a weak population turnover, the disease persists no matter the vaccination coverage and whatever trade-off are assumed in the vaccine (Figure A6a–c).

5. Discussion and Conclusions

When a large proportion of a population becomes immune to a virus, it becomes harder for the disease to spread. This is the core concept underlying the concept of herd immunity [7,8,9]. However, there are numerous individuals who refuse to be vaccinated because of various reasons (health concerns, lack of information, systemic mistrust, see [30]), and some vaccines provide only partial protection from disease or can be only efficient against a few disease variants (see the recent COVID-19 epidemics and the vaccine efficiency and waning of immunity against different variants). Therefore, it is rather common that pathogens face a heterogeneous population of vaccinated and unvaccinated hosts [30], and this has consequences for the evolution of the disease itself [3,15,19]. In this study, we used mathematical modelling approaches (analysis and numerical simulations) to assess the potential population-level impact of using different types of imperfect vaccines to control the burden of a disease in a community. In the first part, we provide a theoretical analysis of the model, including the basic reproduction number ${\mathcal{R}}_0$ and conditions for the stability of the equilibria. We derive the condition to be satisfied regarding the proportion of vaccinated individuals at a steady state in order to attain herd immunity. We express this condition as the critical coverage to be achieved for ${\mathcal{R}}_0<1$.

When the vaccine is developed to prevent infection and stop transmission, our results show that it is possible to eliminate the disease with a strong population turnover if the vaccination coverage is greater than $69.6$% ($48.9$%) with a weak (strong) efficiency of the vaccine. However, when population turnover is weak, we observe damped oscillations, and eradication is possible with a vaccine with high efficiency and coverage greater than $75.6$%. Otherwise, the disease persists and becomes endemic in the community. We highlight here the effect of population turnover as an important first factor in deciding the effectiveness of vaccination campaigns (as suggested in [10,11,13]). For example, with respect to applications to a human population, the turnover can be considered migration in and out of the community since the birth and death rate are usually small and fairly constant. Our results suggest that for a community with strong migration (strong turnover), we can vaccinate individuals coming in, in order to reduce the basic reproduction number. However, if there is a weak migration (weak turnover) as, for example, when flights and travel are restricted, the vaccination strategy should be improved by undertaking a mass vaccination campaign and using a high-efficiency vaccine. A similar reasoning applies to domesticated animals (livestock), with the migration of (potentially vaccinated) individuals between farms influencing the epidemic.

We then analyse more finely the effect of the type of vaccine and its efficiency on disease dynamics. The vaccine can decrease transmission and/or favour the recovery of infected individuals. Disease eradication is possible if the vaccine decreases transmission by 82%, enhances recovery by at least 25%, and a vaccination coverage of 82% is achieved under a strong population turnover. Under weak turnover, maximum vaccine efficiency and coverage are required. Therefore, there is also an interplay between the strength of population turnover and the efficiency of the vaccine (and the property of the vaccine). Finally, we explore the importance of vaccine design is a trade-off between the vaccine efficiency to stop transmission (infection) and disease recovery is expected. We use three trade-off curves and show that the convex ($\nu ={\epsilon }^2$) function is the most desirable when the efficiency of the vaccine is at least 60% under a strong turnover of population. However, under a weak population turnover, the disease cannot be easily eradicated, no matter the vaccination coverage and the efficiency of a combined vaccine. Furthermore, we notice that a smaller vaccination coverage and/or efficiency is needed when using a vaccine designed with a convex trade-off between the above two properties (decrease transmission and favour recovery) than other vaccines (different trade-offs or no trade-off).

Our model has some limitations and advantages compared to previous work in the literature, as we intend here to study the overall behaviour of our model under different schematic scenarios. First, we use, for illustrative purposes, COVID-19 parameters to exemplify the expected threshold for vaccination coverage for a highly transmissible disease. We thus caution here against building precise recommendations (for COVID-19 vaccination) based on our results. Second, our model does not explicitly account for a continuous vaccination (or large vaccination campaigns) of individuals in a community. Vaccination is linked in our model to population turnover, explaining the appearance of periodic oscillations in disease incidence (the honeymoon periods, [7,10,11,28]). Such periodic epidemics occur and are predicted for COVID-19, and may likely be due to immunity waning of the various vaccines against new variants [9]. Third, we use a frequency-dependent transmission, which allows us to derive analytical results in more depth than some previous models but may underestimate the spread of disease and speed of disease dynamics. Thus, to obtain precise predictions regarding vaccination efficiency and campaigns for a given disease, the ad hoc parameters of our models need to be correctly adjusted.

This model contains some general conclusions that are not only applicable to human populations but also domesticated (livestock) and wild animals or even crops. Domesticated animals also require vaccinations (e.g., [31,32]), and our study draws recommendations on the importance of turnover and migration rates in and out of the population. Our results also suggest that in livestock, the type of vaccine can be adjusted depending on the disease, especially if it is desirable that infected animals recover well rather than attempting to prevent any transmission (e.g., [31,32]). Our results may also be relevant to consider for vaccination campaigns of wild endangered animals [33]. In addition, we also suggest that the principles of the model apply to plant (crop) immunisation. To protect plants against invasion of pathogens or pests, one can use different biotic and synthetic chemicals to induce immunity in the plant [34] or protect plants by spraying fungicides [35]. In a field or among fields, some plants will be more resistant than others for a certain period of time. The spray is equivalent to the vaccination, and is, in that case, decoupled from the population turnover, which is the planting/renewal and harvesting/removal of plants. Plant epidemiology modelling has been used to predict the efficiency of imperfect fungicide treatments on the epidemics and on yield [35,36], with results mirroring our own.

In summary, we show that it is possible to achieve disease control by vaccination in a population with a strong turnover, even if we use a weak imperfect vaccine designed to reduce only transmission. However, higher vaccination coverage and a strong efficiency vaccine are necessary to control the disease under weak population turnover. Moreover, a vaccine with a convex trade-off between the efficiency to reduce transmission and enhance recovery is recommendable, along with high vaccination coverage.