# Impact of imperfect vaccine, vaccine trade-off and population turnover on
infectious disease dynamics

Hetsron L. Nyandjo-Bamen  Jean Marie Ntaganda  Aurélien Tellier  Olivier
Menoukeu-Pamen  Department of Mathematics, School of Science, College of
Science and Technology, University of Rwanda, Rwanda  African Institute for
Mathematical Sciences, Ghana  Population Genetics, Department of Life Science
systems, School of Life Sciences, Technical University of Munich 85354
Freising, Germany  IFAM, Department of Mathematical Sciences, University of
Liverpool, United Kingdom

###### Abstract

Vaccination is essential for the management of infectious diseases, many of
which continue to pose devastating public health and economic challenges
across the world. 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 trade-off between these two
properties. Furthermore, population turnover, that is the rate at which
individuals enter and exit the population, is another key factor determining
the epidemiological dynamics. While these factors have yet been studied
separately, we investigate the interplay between the efficiency and property
of an imperfect vaccine and population turnover. We build a mathematical model
with frequency incidence rate, a recovered compartment, and an heterogeneous
host population with respect to vaccination. We first compute the basic
reproduction number $\mathcal{R}_{0}$ and study the global stability of the
equilibrium points. Using a sensitivity analysis, we then assess the most
influential parameters determining the total number of infected and
$\mathcal{R}_{0}$ over time. We derive analytically and numerically conditions
for the vaccination coverage and efficiency to achieve disease eradication
($\mathcal{R}_{0}<1$) assuming different intensity of the population turnover
(weak and strong), vaccine properties (transmission and/or recovery) and
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 trade-off between the vaccine properties. We conclude that
the coverage target for vaccination campaigns should be evaluated based on the
interplay between these factors.

###### keywords:

Imperfect vaccine; Vaccination coverage; Vaccine trade-off; Population
turnover; Mathematical model; Global stability; Sensibility analysis.

††journal: ….

##  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, rinderpest
in most parts of the world, ultimately saving the lives of millions of humans
and animals. A perfect vaccine would keep vaccinated people from becoming
infected when exposed to the pathogen. An imperfect vaccine, one that does not
prevent vaccinated individuals from becoming infected upon pathogen exposure,
may still be beneficial in various ways [4]. 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 [5], or preventing
transmission by infected hosts [6]. As we have seen recently with the epidemic
of Covid-19, imperfect vaccines can be used to reduce the number of infected
individuals, but also to protect individuals at risk of developing the more
lethal form of the infection, especially when the efficiency of vaccination
may be volatile and decreases due to the appearance of new variants of the
virus [7, 8, 9].

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 [10, 11, 12]. However, a continuous influx
of new susceptible, possibly unvaccinated individuals contributes to the
disease’s long-term persistence in the population [13, 14]. A frequent
introduction of pathogen into a partially immune population with intermediate
levels of population immunity can lead to an epidemic of longer duration
and/or higher total number of infectious individuals than the introduction
into a naive population [14]. This phenomenon is named as ”epidemic
enhancement” [14]. 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 [15]. In human but also
domesticated animals, population turnover takes the form of immigration and
emigration in and out of the population, as well as birth and death of
individuals. The turnover is an often neglected factor in epidemiology when
generalizing predictions of modelling from human to domesticated and animal
populations.

Moreover, a second parameter of importance in studying the efficiency of
vaccination strategies, is the existence of biological trade-offs in
epidemiology. 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 Gandon et al. [6], it is predicted that
vaccines affecting disease transmission may lead to a decrease of parasite
virulence, while other types of vaccines (reducing within-host growth rate)
may lead to an increase of parasite virulence, and thus the counter-effect of
a worst epidemiological outcome. Interestingly, much work has been devoted to
generate precise predictions for virulence evolution in known parasites by
incorporating empirical characterizations of vaccine effects into models
capturing the epidemiological details of a given system [16, 17, 18]. In
contrast, biochemical and immunological trade-offs of the vaccine itself have
received little attention. We 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 generalize the definition of imperfect vaccines as
providing partial protection (non-maximal efficiency) against infection
(decreasing transmission), partially enhancing (not fully) recovery of
infected individuals, and possible trade-off between these two properties.
There has been remarkably little work done to generally assess how the
interplay between different vaccine properties, trade-offs, and vaccination
strategies influences the burden of the epidemic in an heterogeneous community
with imperfect vaccination.

The aim of this study is therefore to assess, through mathematical modelling,
whether the use of vaccines that decrease the infection is more efficient to
eradicate the disease in an 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 generalize
our results to animal populations.
The paper is organized 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 a local and global stability of the equilibrium points of
the model (Section 3). We then provide a numerical sensitivity of the model
and examples of numerical analyses for different parameter values describing
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 but also
domesticated animal species for which turnover rates represent different end
of a continuum.

##  2 Model formulation

The formulation of the model is based on compartmental modeling [19], which
consists in 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 dynamic of an infectious disease during and after the vaccination
campaign, we modify the model formulated in [6] by adding a recovered
compartment and we consider a frequency-dependent disease transmission
(incidence rate). The model takes in to account only host-to-host transmission
of the disease. Since many vaccines do not guaranty a perfect immunity, we
consider an heterogeneous host community 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. And among the partially
resistant hosts, we find the uninfected ($S_{2}$) and the 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 [20]. Thus, the total population at time $t$, $N(t)$ is
given by

$N(t)=S_{1}(t)+S_{2}(t)+I_{1}(t)+I_{2}(t)+R(t).$

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 constant rate, $\theta$. Among these new uninfected, a
proportion, $p$, is partially immune due to the vaccination, while the
remaining proportion $1-p$ is completely vulnerable to the parasite.
Uninfected, infected and recovered hosts die naturally at a rate, $\mu$ and
infected hosts suffer additional mortality due to the virulence of the
parasite. Since host resistance may reduce the impact of parasite [6], we
assume the virulence of the parasite on fully susceptible hosts, $d_{1}$, is
greater than the one on partially resistant hosts, $d_{2}.$ Uninfected hosts
become infected with the forces of infection
$\lambda_{1}(t)=\beta_{11}\dfrac{I_{1}(t)}{N(t)}+\beta_{12}\dfrac{I_{2}(t)}{N(t)}$
and
$\lambda_{2}(t)=\beta_{21}\dfrac{I_{1}(t)}{N(t)}+\beta_{22}\dfrac{I_{2}(t)}{N(t)}$
when they are fully susceptible or partially resistant, respectively. And
since the resistance can decrease the probability of becoming infected [6], we
assume $\beta_{21}\leqslant\beta_{11}$ and $\beta_{22}\leqslant\beta_{12}$.
Recovery rates may differ between the fully susceptible, $\gamma_{1}$ and the
partially resistant host, $\gamma_{2}$. The schematic diagram of the model is
as shown in Figure 1.

Figure 1: Schematic diagram of the epidemiological model with imperfect
vaccination.

Mathematically, the model is as follows:

$\left\\{\begin{array}[]{llll}\dfrac{\mathrm{d}S_{1}}{\mathrm{d}t}&=&\theta(1-p)-\lambda_{1}(t)S_{1}(t)-\mu
S_{1}(t),\\\ \\\ \dfrac{\mathrm{d}S_{2}}{\mathrm{d}t}&=&\theta
p-\lambda_{2}(t)S_{2}(t)-\mu S_{2}(t),\\\ \\\
\dfrac{\mathrm{d}I_{1}}{\mathrm{d}t}&=&\lambda_{1}(t)S_{1}(t)-(\mu+\gamma_{1}+d_{1})I_{1}(t),\\\
\\\
\dfrac{\mathrm{d}I_{2}}{\mathrm{d}t}&=&\lambda_{2}(t)S_{2}(t)-(\mu+\gamma_{2}+d_{2})I_{2}(t),\\\
\\\
\dfrac{\mathrm{d}R}{\mathrm{d}t}&=&\gamma_{1}I_{1}(t)+\gamma_{2}I_{2}(t)-\mu
R(t).\end{array}\right.$ (1)

A summary of the biological significance of the model’s parameters (1) is
given in Table 1.

Table 1: Description and value of the model’s 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 form $I_{1}$ to $S_{1}$ | $day^{-1}$ | variable | Assumed
$\beta_{12}$ | Transmission rate form $I_{1}$ to $S_{2}$ | $day^{-1}$ | variable | Assumed
$\beta_{21}$ | Transmission rate form $I_{2}$ to $S_{1}$ | $day^{-1}$ | variable | Assumed
$\beta_{22}$ | Transmission rate form $I_{2}$ to $S_{2}$ | $day^{-1}$ | variable | Assumed
$d_{1}$ | Mortality rate due to infection of $S_{1}$ | $day^{-1}$ | $0.0008$ | [12]
$d_{2}$ | Mortality rate due to infection of $S_{2}$ | $day^{-1}$ | $0.0001$ | [12]
$\gamma_{1}$ | Recovery rate of $I_{1}$ | $day^{-1}$ | $0.1$ | [12]
$\gamma_{2}$ | Recovery rate of $I_{2}$ | $day^{-1}$ | $0.13$ | [12]

##  3 Mathematical analysis

###  3.1 Basic properties

First, we study the basic characteristics of the system solutions: the
existence, positivity 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 Positivity of solutions

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

######  Theorem 3.1.

For any initial condition
$(t_{0}=0,X_{0}=(S_{1}(0),S_{2}(0),I_{1}(0),I_{2}(0),R(0))\in\mathbb{R}^{5}_{+})$,
and for $T\in]0,+\infty]$, the maximal solution
$([0,T[,X=(S_{1}(t),S_{2}(t),I_{1}(t),I_{2}(t),R(t)))$ of the Cauchy problem
associated to the system (1) is non-negative.

###### Proof.

Let

$\Delta=\\{\tilde{t}\in[0,T[\hskip 5.69046pt|\hskip
5.69046ptS_{1}(t)>0,S_{2}(t)>0,I_{1}(t)>0,I_{2}(t)>0\hskip
5.69046pt\text{and}\hskip 5.69046ptR(t)>0\hskip 5.69046pt\forall
t\in[0,\tilde{t}[\\}.$

By continuity of the functions $S_{1}$, $S_{2}$, $I_{1}$, $I_{2}$ and $R$, one
can see that $\Delta\neq\emptyset$. Let $\tilde{T}=\text{sup}\Delta$ be the
supremum of $\Delta$. Now, we want to prove that $\tilde{T}=T$.

Assume $\tilde{T}<T,$ then we have that $S_{1}$, $S_{2}$, $I_{1}$, $I_{2}$ and
$R$ are simultaneously positive on $]0,\tilde{T}[$. Then, at least one of the
following conditions is satisfied at time $\tilde{T}$: $S_{1}(\tilde{T})=0$
and $\dfrac{\mathrm{dS_{1}(\tilde{T})}}{\mathrm{d}t}\leqslant 0$,
$S_{2}(\tilde{T})=0$ and
$\dfrac{\mathrm{dS_{2}(\tilde{T})}}{\mathrm{d}t}\leqslant 0$,
$I_{1}(\tilde{T})=0$ and
$\dfrac{\mathrm{dI_{1}(\tilde{T})}}{\mathrm{d}t}\leqslant 0$,
$I_{2}(\tilde{T})=0$ and
$\dfrac{\mathrm{dI_{2}(\tilde{T})}}{\mathrm{d}t}\leqslant$ or $R(\tilde{T})=0$
and $\dfrac{\mathrm{dR(\tilde{T})}}{\mathrm{d}t}\leqslant 0$.

Suppose that $S_{1}(\tilde{T})=0$ and
$\dfrac{\mathrm{dS_{1}(\tilde{T})}}{\mathrm{d}t}\leqslant 0$, then we deduce
from the first equation of system (1) that

$\displaystyle\dfrac{\mathrm{d}S_{1}(\tilde{T})}{\mathrm{d}t}$
$\displaystyle=\theta(1-p)-\lambda_{1}(\tilde{T})S_{1}(\tilde{T})-\mu
S_{1}(\tilde{T})$ $\displaystyle=\theta(1-p)>0$

which is a contradiction to the previous claim that
$\dfrac{\mathrm{dS_{1}(\tilde{T})}}{\mathrm{d}t}\leqslant 0$.

We can use the similar argument for all the remaining state variables. Then,
$\tilde{T}=T$ and consequently the maximal solution
$(S_{1}(t),S_{2}(t),I_{1}(t),I_{2}(t),R(t))$ of the Cauchy problem related
with system (1) is positive. ∎

Therefore, the variables of the system (1) are positive for all time $t>0$. In
other terms, solutions of the system (1) with non-negative initial conditions
will stay positive for all $t>0$.

We now prove a useful lemma.

####  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 realistic that the
total number of individuals does not tend towards infinity. The following
result yields:

######  Lemma 3.2.

The closed set

$\Omega=\bigg{\\{}(S_{1}(t),S_{2}(t),I_{1}(t),I_{2}(t),R(t))\in\mathbb{R}^{5}_{+},\hskip
5.69046ptN(t)\leqslant\dfrac{\theta}{\mu}\bigg{\\}}$

is positively invariant and attracting for the system (1).

###### Proof.

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

$\dfrac{\mathrm{d}N}{\mathrm{d}t}=\theta-\mu
N-d_{1}I_{1}-d_{2}I_{2}\leqslant\theta-\mu N.$

Integrating both sides of the expression above, we deduce that

$N(t)\leqslant\dfrac{\theta}{\mu}+\bigg{(}N(0)-\dfrac{\theta}{\mu}\bigg{)}e^{-\mu
t},\quad\forall t\geqslant 0,$ (2)

where $N(0)$ is the value of $N(t)$ at the beginning.

We deduce that if $N(0)\leqslant\dfrac{\theta}{\mu}$, then $0\leqslant
N(t)\leqslant\dfrac{\theta}{\mu}$, $\forall t\geqslant 0$ and $\Omega$ is
positively invariant. If $N(0)\geqslant\dfrac{\theta}{\mu}$, then from (2) the
total population decreases and the solutions enters $\Omega$. Hence $N(t)$ is
bounded as $t$ →$\infty$, which means that $\Omega$ is attracting. ∎

######  Remark 3.1.

In finite dimension every maximal solution of a Cauchy problem is global in a
compact set. Then, every maximal solution of the model system (1) is global.

Therefore the solutions of our model are considered epidemiologically and
mathematically well posed in $\Omega$ .

###  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 a state in the population without any infection.
The disease-free equilibrium is deduced from the resolution of the system of
equations in (1) by taking $I_{1}=0$ and $I_{2}=0$. Thus, the disease-free
equilibrium for model (1) satisfies the following system of equations:

$\left\\{\begin{array}[]{ll}\theta(1-p)-\mu S_{1}^{0}=0,\\\ \theta p-\mu
S_{2}^{0}=0.\end{array}\right.$ (3)

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

$Q^{0}=(S_{1}^{0},S_{2}^{0},0,0,0),$

where $S_{1}^{0}=\dfrac{\theta(1-p)}{\mu}$, $S_{2}^{0}=\dfrac{\theta p}{\mu}$
and $N^{0}=S_{1}^{0}+S_{2}^{0}=\dfrac{\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 infectivity period when it is
introduced into a population of susceptible individuals. We study the
stability of the equilibrium through the next generation operator [21, 22].
Recalling the notations in [22] for model (1), the matrices $\mathcal{F}$ of
the new infection and $\mathcal{V}$ of the remaining transfer terms at the DFE
for are given by

$\mathcal{F}=\begin{bmatrix}\beta_{11}\dfrac{S_{1}I_{1}}{N}+\beta_{12}\dfrac{S_{1}I_{2}}{N}\\\
\\\
\beta_{21}\dfrac{S_{2}I_{1}}{N}+\beta_{22}\dfrac{S_{2}I_{2}}{N}\end{bmatrix}\hskip
5.69046pt\text{and}\hskip
5.69046pt\mathcal{V}=\begin{bmatrix}(\mu+\gamma_{1}+d_{1})I_{1}\\\ \\\
(\mu+\gamma_{2}+d_{2})I_{2}\end{bmatrix}.$

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

$F=\begin{bmatrix}\beta_{11}\dfrac{S_{1}^{0}}{N^{0}}&\beta_{12}\dfrac{S_{1}^{0}}{N^{0}}\\\
\\\
\beta_{21}\dfrac{S_{2}^{0}}{N^{0}}&\beta_{22}\dfrac{S_{2}^{0}}{N^{0}}\end{bmatrix}\hskip
5.69046pt\text{and}\hskip 5.69046ptV=\begin{bmatrix}\mu+\gamma_{1}+d_{1}&0\\\
\\\ 0&\mu+\gamma_{2}+d_{2}\end{bmatrix}.$ (4)

Then,

$FV^{-1}=\begin{bmatrix}\dfrac{\beta_{11}S_{1}^{0}}{N^{0}(\mu+\gamma_{1}+d_{1})}&\dfrac{\beta_{12}S_{1}^{0}}{N^{0}(\mu+\gamma_{2}+d_{2})}\\\
\\\
\dfrac{\beta_{21}S_{2}^{0}}{N^{0}(\mu+\gamma_{1}+d_{1})}&\dfrac{\beta_{22}S_{2}^{0}}{N^{0}(\mu+\gamma_{2}+d_{2})}\end{bmatrix},$

and the reproduction number of model system (1) is

$\displaystyle\mathcal{R}_{0}=\rho(FV^{-1})=$
$\displaystyle\dfrac{1}{2}\Big{[}\dfrac{S_{1}^{0}}{N^{0}}\mathcal{R}_{0,11}+\dfrac{S_{2}^{0}}{N^{0}}\mathcal{R}_{0,22}+\sqrt{\Big{(}\dfrac{S_{1}^{0}}{N^{0}}\mathcal{R}_{0,11}-\dfrac{S_{2}^{0}}{N^{0}}\mathcal{R}_{0,22}\Big{)}^{2}+4\dfrac{S_{1}^{0}}{N^{0}}\dfrac{S_{2}^{0}}{N^{0}}\mathcal{R}_{0,12}\mathcal{R}_{0,21}}\Big{]},$
$\displaystyle\mathcal{R}_{0}=$
$\displaystyle\dfrac{1}{2}\Big{[}(1-p)\mathcal{R}_{0,11}+p\mathcal{R}_{0,22}+\sqrt{\Big{(}(1-p)\mathcal{R}_{0,11}-p\mathcal{R}_{0,22}\Big{)}^{2}+4p(1-p)\mathcal{R}_{0,12}\mathcal{R}_{0,21}}\Big{]},$
(5)

where $\dfrac{S_{1}^{0}}{N^{0}}=1-p$ (respectively
$\dfrac{S_{2}^{0}}{N^{0}}=p$) is the proportion of susceptible individuals
that have not been vaccinated (respectively have been vaccinated) at the DFE
$Q^{0}$. Similarly, we define
$\mathcal{R}_{0,11}=\dfrac{\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 the interaction with the unvaccinated
population. Also $\mathcal{R}_{0,12}=\dfrac{\beta_{12}}{\mu+\gamma_{1}+d_{1}}$
represents the average number of secondary cases generated by a vaccinated
infected in the unvaccinated part of the population,
$\mathcal{R}_{0,21}=\dfrac{\beta_{21}}{\mu+\gamma_{2}+d_{2}}$ is the average
number of secondary cases generated by an unvaccinated infected in the
vaccinated part of the population, and
$\mathcal{R}_{0,22}=\dfrac{\beta_{22}}{\mu+\gamma_{2}+d_{2}}$ represents the
average number of secondary cases generated by a vaccinated infected in the
vaccinated part of the population. Further, $\rho(FV^{-1})$ is the spectral
radius of $FV^{-1}$.

######  Remark 3.2.

From the expression of the reproduction number $\mathcal{R}_{0}$ in (5), we
deduce that
$\mathcal{R}_{0}\geq(1-p)\mathcal{R}_{0,11}\vee p\mathcal{R}_{0,22}$. Moreover
using (5) for $p=0$ (all new hosts are not vaccinated),
$\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 [22].

######  Lemma 3.3.

The DFE $Q^{0}$ of the system (1) is locally asymptotically stable whenever
$\mathcal{R}_{0}<1$ and unstable whenever $\mathcal{R}_{0}>1$.

The biological meaning of Lemma 3.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. Global asymptotic stability (GAS) of
the DFE is required to better control the disease. In addition, 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 [23].

######  Theorem 3.4.

If $\mathcal{R}_{0}\leqslant 1$, the DFE $Q^{0}$ of the system (1) is GAS in
$\Omega$. If $\mathcal{R}_{0}>1$, $Q^{0}$ is unstable, the system (1) is
uniformly persistent and there exists at least one endemic equilibrium in the
interior of $\Omega$.

###### Proof.

See A. ∎

As a consequence of the meaning of Theorem 3.4 and Remark 3.2, we can
confidently deduce that the disease can be eradicated from the host community
if the value of $\mathcal{R}_{0}$ can be reduced to less than the unity,
independently of whether individuals introduced in 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}[]{llll}\theta(1-p)-\beta_{11}\dfrac{S_{1}^{*}I_{1}^{*}}{N^{*}}-\beta_{12}\dfrac{S_{1}^{*}I_{2}^{*}}{N^{*}}-\mu
S_{1}^{*}=0,\\\ \\\ \theta
p-\beta_{21}\dfrac{S_{2}^{*}I_{1}^{*}}{N^{*}}-\beta_{22}\dfrac{S_{2}^{*}I_{2}^{*}}{N^{*}}-\mu
S_{2}^{*}=0,\\\ \\\
\beta_{11}\dfrac{S_{1}^{*}I_{1}^{*}}{N^{*}}+\beta_{12}\dfrac{S_{1}^{*}I_{2}^{*}}{N^{*}}-(\mu+\gamma_{1}+d_{1})I_{1}^{*}=0,\\\
\\\
\beta_{21}\dfrac{S_{2}^{*}I_{1}^{*}}{N^{*}}+\beta_{22}\dfrac{S_{2}^{*}I_{2}^{*}}{N^{*}}-(\mu+\gamma_{2}+d_{2})I_{2}^{*}=0,\\\
\\\ \gamma_{1}I_{1}^{*}+\gamma_{2}I_{2}^{*}-\mu R^{*}=0.\end{array}\right.$
(6)

Given the complexity of the system (6), 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 its size is large. However, to prove the existence of $Q^{*}$, we can
rewrite the system (6) as a fixed point problem and use Theorem 2.1 in [24].
To do this, we solve the system (6). After algebraic manipulations of this
system, we obtain:

$R^{*}=\dfrac{\gamma_{1}I_{1}^{*}+\gamma_{2}I_{2}^{*}}{\mu}$,
$S_{1}^{*}=\dfrac{\theta(1-p)N^{*}}{\beta_{11}I_{1}^{*}+\beta_{12}I_{2}^{*}+\mu
N^{*}}$, $S_{2}^{*}=\dfrac{\theta
pN^{*}}{\beta_{21}I_{1}^{*}+\beta_{22}I_{2}^{*}-d_{1}I_{1}^{*}-d_{2}I_{2}^{*}+\theta}$,

$I_{1}^{*}=\dfrac{\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}(I^{*})$
and

$I_{2}^{*}=\dfrac{\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}(I^{*})$
with $I^{*}=(I_{1}^{*},I_{2}^{*}).$

Then, the endemic equilibrium are the fixed points of $H$ given by $I=H(I)$
where $I=(I_{1},I_{2})$. By definition, $H$ is continuous, monotonously non
decreasing and strictly sublinear. $H$ is also a bounded function which maps
the non negative orthant $\Omega$ into itself. Morever, $H(0)=0$ by definition
and the jacobian of $H$ at the zero, $H^{{}^{\prime}}(0)$, exists and is
irreducible since

$H^{{}^{\prime}}(0)=\begin{bmatrix}\beta_{11}a_{1}&\beta_{12}a_{1}\\\ \\\
\beta_{21}a_{2}&\beta_{22}a_{2}\end{bmatrix}=FV^{-1},$

where $a_{1}=\dfrac{1-p}{\mu+\gamma_{1}+d_{1}}$ and
$a_{2}=\dfrac{p}{\mu+\gamma_{2}+d_{2}}.$

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

######  Proposition 3.1.

The system (1) has only one endemic equilibrium whenever $\mathcal{R}_{0}>1$.

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

######  Theorem 3.5.

If $\mathcal{R}_{0}>1$, the endemic equilibrium $Q^{*}$ is GAS in $\Omega$.

###### Proof.

See B. ∎

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
can occur with some diseases 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
who lack 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 the spread of disease. To compute the herd
immunity threshold associated with the model (1), we set the reproduction
number, $\mathcal{R}_{0}$ to one and solve for $p=\dfrac{S_{2}^{0}}{N^{0}}$
which is the proportion of susceptible individuals which have been vaccinated
at the DFE, $Q^{0}.$ Then we have,

$\displaystyle\mathcal{R}_{0}=1\Longleftrightarrow$
$\displaystyle\big{[}2-\mathcal{R}_{0,11}+(\mathcal{R}_{0,11}-\mathcal{R}_{0,22})p\big{]}^{2}=\big{[}\mathcal{R}_{0,11}-(\mathcal{R}_{0,11}+\mathcal{R}_{0,11})p\big{]}^{2}+4p(1-p)\mathcal{R}_{0,12}\mathcal{R}_{0,21}$
$\displaystyle\Longleftrightarrow$
$\displaystyle\big{[}(\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}\big{]}p^{2}+\big{[}2(2-\mathcal{R}_{0,11})(\mathcal{R}_{0,11}-\mathcal{R}_{0,22})$
$\displaystyle+2\mathcal{R}_{0,11}(\mathcal{R}_{0,11}+\mathcal{R}_{0,22})-4\mathcal{R}_{0,12}\mathcal{R}_{0,21}\big{]}p+(2-\mathcal{R}_{0,11})^{2}-\mathcal{R}_{0,11}^{2}=0.$

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

$Q(p)=Ap^{2}+Bp+C,$ (7)

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, the conditions
for $Q(p)$ to have positive real roots are determined below. For this purpose,
we now perform a case analysis to determine the positive real zeros of $Q$.

Let $\Delta=B^{2}-4AC$ be the discriminant of the equation $Q(p)=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\neq 0.$

Case 2

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

$p_{c_{0}}=-\dfrac{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\neq 0$ and $\Delta>0$. Then

$p_{c_{1}}=\dfrac{-B-\sqrt{\Delta}}{2A}\hskip 5.69046pt\text{and}\hskip
5.69046ptp_{c_{2}}=\dfrac{-B+\sqrt{\Delta}}{2A}$

are the real roots of $Q$.

Moreover, if $A>0$, then

$\left\\{\begin{array}[]{ll}p_{c_{1}}>0\hskip 5.69046pt\text{if and only
if}\hskip 5.69046pt\sqrt{\Delta}<-B,\\\ p_{c_{2}}>0\hskip 5.69046pt\text{if
and only if}\hskip 5.69046pt\sqrt{\Delta}>B.\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}[]{ll}p_{c_{1}}>0\hskip 5.69046pt\text{if and only
if}\hskip 5.69046pt\sqrt{\Delta}>-B,\\\ p_{c_{2}}>0\hskip 5.69046pt\text{if
and only if}\hskip 5.69046pt\sqrt{\Delta}<B.\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$).

Theorem 3.4 and Theorem 3.5 can be combined to give the following result:

######  Corollary 3.1.

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

The implication of Corolloary 3.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 individuals vaccinated satisfies one of these conditions:

  1. 1.

$p>p_{c}$, if $A=0$, $B>0$ and $C<0$;

  2. 2.

$p\in[0,p_{c}[$, if $A=0$, $B>0$ and $C>0$;

  3. 3.

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

  4. 4.

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

  5. 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. 6.

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

  7. 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. 1.

$p\in[0,p_{c}[$, if $A=0$, $B>0$ and $C<0$;

  2. 2.

$p>p_{c}$, if $A=0$, $B>0$ and $C>0$;

  3. 3.

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

  4. 4.

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

  5. 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. 6.

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

  7. 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$).

To conclude on the analytical part, the eradication of a disease is
conditioned by the proportion of vaccinated individuals, this vaccination
coverage threshold is called the critical vaccination proportion ($p_{c}$). In
some cases, there is one critical proportion which determines whether the
basic reproduction number,$\mathcal{R}_{0}$, is less than one or not. In other
cases, two critical proportions are found and which define three different
dynamics: disease eradication when $\mathcal{R}_{0}<1$, endemic disease
dynamics when $\mathcal{R}_{0}>1$ with presence or absence of epidemiological
oscillations in the number of infected. In the latter case of two thresholds,
the analytical results derive above do not allow the prediction of occurrence
of the 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 behavior of model (1), we use parameter values for the
mortality rates, $d_{1}$, $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 [12]), and vary
the values of other parameters as described in Table 1.

###  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 [25].
The sensitivity analysis of the model parameters is carried out to measure the
correlation between the model’s parameters (1) and 1) the total number of
infected individuals ($I_{1}+I_{2}$), and 2) the threshold parameter
$\mathcal{R}_{0}$. The analysis is performed by using the Latin Hypercube
Sampling (LHS) technique and partial rank correlation coefficients (PRCCs)
[25]. In our analysis, 1,000 model simulations are performed by running the
model for 200 time steps (equivalent to 200 days) and number of infected are
recorded at time points 50, 100 and 200. To perform the sensitivity analysis,
each parameter has a parameter range defines by the maximum (respectively the
minimum) being $50\%$ greater (respectively less) than its baseline (values in
Table 3, 4, 5, 6). We then divide each parameter range into 1,000 equally
large sub-intervals, and draw a value per parameter within that interval using
a Uniform draw. By this mean we obtain a uniform distribution of 1,000
parameter values for each parameter. The parameter space (or LHS matrix) has
dimension of length 11 with each dimension specifying an uncertain parameter
vector of length 1,000. The base parameter values are chosen to define several
scenarios of interest regarding the intensity of the turnover (weak and
strong) and efficiency of the vaccine (weak and strong). In PRCC analysis, the
parameters with the larger positive or negative PRCC values ($>0.5$ or
$<-0.5$) and with corresponding 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 Tables 3, 4, 5, 6
in C .

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

| Scenarios | Total Infected: $I_{1}+I_{2}$
---|---
|  $t=50$ days |  $t=100$ days |  $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}(+)$

Based on the results from Tables 3, 4, 5, 6, we provide in table 2, a summary
of the the parameters that significantly affect the number of infected.
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. This suggests that an effective
control strategy should aim to limit significantly the immigration of new
hosts in the population (to decrease $\theta$) and improve the treatment of
infected people (to increase $\gamma_{1}$). We then proceed to a similar
analysis with $\mathcal{R}_{0}$, and summarize 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).

(a) Strong turnover and weak efficiency

(b) Strong turnover and strong efficiency

(c) Weak turnover and weak efficiency

(d) Weak turnover and strong efficiency

Figure 2: PRCCs describing the impact of model’s parameters on
$\mathcal{R}_{0}$ of the model (1) with respect to some scenarios. The range
of the parameters in (a) (respectively in (b), (c) and (d)) is the same as
given on Table 3(respectively on Table 4,5,6).

###  4.2 Interplay between vaccine efficiency and population turnover

We now study the effect of population turn-over 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 3.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\%$).

(a)

(b)

(c)

(d)

Figure 3: Epidemiological dynamics with the initial conditions
$S_{1}(0)=1000$, $S_{2}(0)=700$, $I_{1}(0)=200$, $I_{2}(0)=80$, $R(0)=20$ for
various scenarios assuming the parameters $\beta_{11}=0.35$,
$\beta_{12}=0.28$, $p=0.5$ and strong population turnover ($\theta=1000$,
$\mu=0.09$). We present under weak vaccine efficiency
($\beta_{21}=0.175,\beta_{22}=0.14$), the number of (a) uninfected and (b)
infected individuals. We present under strong vaccine efficiency
($\beta_{21}=0.035,\beta_{22}=0.028$) the number of (c) uninfected and (d)
infected individuals. Others parameters values are as in Table 1.

####  4.2.1 Strong population turnover

The epidemiological dynamics in Figure3(b) under strong turnover and weak
vaccine efficiency ($\mathcal{R}_{0}=1.2352$) shows that the dynamics reaches
the endemic disease equilibrium. Furthermore if $p$ takes 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 3.1). So to eradicate the disease under strong population
turnover and weak efficiency of the vaccine, a minimum vaccination rate is
needed and defined by $p_{1}$. Under strong turnover and strong efficiency
(Figure3(d), with $\mathcal{R}_{0}=0.9808$) the disease becomes extinct.
Furthermore if the parameter $p$ 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$. So 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.

(a)

(b)

(c)

(d)

Figure 4: Simulation of model (1) at the initial conditions $S_{1}(0)=1000$,
$S_{2}(0)=700$, $I_{1}(0)=200$, $I_{2}(0)=80$, $R(0)=20$ when $\theta=10$,
$\beta_{11}=0.35$, $\beta_{12}=0.28$, $\beta_{21}=0.175,\beta_{22}=0.14$,
$\mu=0.0009$, $p=0.5$, (a) Uninfected individuals in weak turnover and weak
efficiency scenario and (b) Infected individuals in weak turnover and weak
efficiency scenario. When $\theta=1000$, $\beta_{11}=0.35$, $\beta_{12}=0.28$,
$\beta_{21}=0.035,\beta_{22}=0.028$, $\mu=0.0009$, $p=0.5$,(c) Uninfected
individuals in weak turnover and strong efficiency scenario and (d) Infected
individuals in weak turnover and strong efficiency scenario. Others parameters
values are as in Table 1.

####  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 exhibits damped oscillations (recurring outbreaks) before stabilizing
at the endemic state with disease persistence (Figure4(b) with
$\mathcal{R}_{0}=2.2551$, Figure4(d) 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 are available from new recruitment into the population and
recovered individuals loosing their immunity (so-called waning immunity). This
phenomenon was also described in [11, 14, 26, 27], and the effect of turnover
and waning immunity is specifically described in [11, 14].
With respect to the control of the disease, under weak vaccine efficiency, $p$
can take any value between $0$ and $1$, the basic reproduction number is
always greater than $1$ (Figure4(b) with $\mathcal{R}_{0}=2.2551$). In
contrast, when vaccine efficiency is strong, three cases occur Figure4(d)
(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 values between
$p_{3}$ and $p_{4}$ (with $p_{4}\approx 0.756$), the basic reproduction
number, $\mathcal{R}_{0}$, is greater than $1$ but there are no periodic
oscillations. And for $p\in[p_{4},1]$, the basic reproduction number,
$\mathcal{R}_{0}$, is less than $1$, and disease becomes extinct. Note that
between $p_{3}$ and $p_{4}$, the behavior 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,
a high vaccination coverage (more than $75.6\%$ of the new host individuals)
is needed. Our results extend those in [28] showing that it is feasible to
control 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 (eq. 5). The simulation results provide examples of the analytical
expressions obtained in eq. 7 following the Corollary 3.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:

$\displaystyle\beta_{21}=(1-\varepsilon)\beta_{11},\beta_{22}=(1-\varepsilon)\beta_{12},\hskip
2.84544pt\text{and}\hskip 2.84544pt\gamma_{1}=(1-\nu)\gamma_{2},$ (8)

where $0\leq\varepsilon\leq 1$ represents the effect of the vaccine on the
transmission and $0\leq\nu\leq 1$ represents the effect of the vaccine on the
ability of being recovered. Substituting the rescaled expressions for (8) into
the model (1), one deduces that the basic reproduction number the model (1)
can be rewritten as:

$\displaystyle\mathcal{R}_{0}=$
$\displaystyle\dfrac{1}{2}\Big{[}(1-p)\mathcal{R}_{0,11}+p\mathcal{R}_{0,22}+\sqrt{\Big{(}(1-p)\mathcal{R}_{0,11}-p\mathcal{R}_{0,22}\Big{)}^{2}+4p(1-p)\mathcal{R}_{0,12}\mathcal{R}_{0,21}}\Big{]},$
(9)

with $\mathcal{R}_{0,11}=\dfrac{\beta_{11}}{\mu+(1-\nu)\gamma_{2}+d_{1}}$,
$\mathcal{R}_{0,12}=\dfrac{\beta_{12}}{\mu+(1-\nu)\gamma_{2}+d_{1}}$,
$\mathcal{R}_{0,21}=\dfrac{(1-\varepsilon)\beta_{11}}{\mu+\gamma_{2}+d_{2}}$
and
$\mathcal{R}_{0,22}=\dfrac{(1-\varepsilon)\beta_{12}}{\mu+\gamma_{2}+d_{2}}$.
Simulations are carried out to assess the interplay of the type of vaccine and
the population turnover.

(a)

(b)

(c)

Figure 5: Contour plots of the basic reproduction number ($\mathcal{R}_{0}$)
of the model (1) with a strong population turnover as a function of (a)
vaccination coverage, $p$, and vaccine efficiency on disease transmission,
$\varepsilon$ (with fixed $\nu=0.5$); (b) vaccination coverage, $p$, and
vaccine efficiency on recovery, $\nu$ (with fixed $\varepsilon=0.5$); and (c)
vaccine efficiency on recovery, $\nu$, and vaccine efficiency on transmission,
$\varepsilon$ (with fixed $p=0.5$). The parameters are $\theta=1000$,
$\beta_{11}=0.35$, $\beta_{12}=0.28$, $\beta_{21}=0.175,\beta_{22}=0.14$,
$\mu=0.09$, $d_{1}=0.0008$, $d_{2}=0.0001$, $\gamma_{1}=0.065$,
$\gamma_{2}=0.13$.

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 5(a)), and if the vaccine is designed to only decrease the
transmission by $80\%$ (i.e. $\varepsilon=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 5(a)). On the other hand,
the value of the reproduction number decreases as coverage increases and
efficiency of the vaccine favoring recovery decreases (Figure 5(b)). 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 5(b)). In Figure 5(c), 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 7, in which
$\mathcal{R}_{0}$ is a function of $\varepsilon$, $\nu$ and $p$. The use of a
vaccine with a combined efficiency (decreasing transmission and favouring
recovery) can be associated to the vaccination coverage in order to achieve
the elimination of the disease. For example, with a vaccination coverage of
$20\%$ ($p=0.2$), it is not possible to eliminate the disease no matter the
combined efficiency of the vaccine (Figure 8), 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
8).

The above results change dramatically 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 9(a)), but a
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 to enhance recovery,
no matter the vaccination coverage (Figure 9(b)), or 2) if the efficiency of
the vaccine is combined but vaccination coverage is $p=0.5$ (Figure 10). The
general results of $\mathcal{R}_{0}$ as a function of $\varepsilon$, $\nu$ and
$p$ demonstrate that under weak population turnover, disease eradication
requires a very strong efficiency of the vaccine and a high coverage (Figure
11).

###  4.4 Interplay between vaccine efficiency trade-off and population
turnover

So 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 the transmission and on the recovery. We
assume three possible trade-off curves: convex($\nu=\varepsilon^{2}$),
concave($\nu=\sqrt{\varepsilon}$) or linear($\nu=\varepsilon$). Under a strong
population turnover, assuming a vaccine of at least $60\%$ of efficiency,
disease eradication can be achieved ($\mathcal{R}_{0}<1$) if the coverage is
at least $65\%$ under a convex trade-off (Figure 6(a)), at least $80\%$ under
a concave trade-off (Figure 6(b)) and at least $75\%$ under a linear trade-off
(Figure 6(c)). Imposing vaccine trade-off affects therefore the shape of the
$\mathcal{R}_{0}$ curves in Figure 6(a), 6(b), 6(c) compared to Figures 5(a)
and 5(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 (Figures 12(a), 12(b) and 12(c)).

(a)

(b)

(c)

Figure 6: Contour plots of the basic reproduction number ($\mathcal{R}_{0}$)
of the model (1) with a strong population turnover as a function of vaccine
coverage, $p$, and vaccine efficiency on the transmission, $\varepsilon$ when:
(a) $\nu=\varepsilon^{2}$ (convex relationship); (b) $\nu=\sqrt{\varepsilon}$
(concave relationship); (c)$\nu=\varepsilon$ (linear relationship). The
parameters are $\theta=1000$, $\beta_{11}=0.35$, $\beta_{12}=0.28$,
$\beta_{21}=0.175,\beta_{22}=0.14$, $\mu=0.09$, $d_{1}=0.0008$,
$d_{2}=0.0001$, $\gamma_{1}=0.065$, $\gamma_{2}=0.13$.

##  5 Discussion and Conclusion

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 [10, 11, 12]. However, there are numerous individuals
who refuse to be vaccinated because of various reasons (health concerns, lack
of information, systemic mistrust, …), and some vaccines provide only partial
protection from disease or can be only efficient against 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 an heterogeneous population of vaccinated and unvaccinated
hosts, and this has consequences for the evolution of the disease itself [6,
17, 20]. In this study, we used mathematical modelling approaches (analysis
and numerical simulations) to assess the potential population-level impact of
the using different types of imperfect vaccines to control the burden of a
disease in a community. In a 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 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
result based on Covid-19 parameter estimates show that it is possible to
eliminate the disease with a strong population turnover if the vaccination
coverage is greater than $69.6$% (respectively $48.9$%) with a weak
(respectively 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 a 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 [13,
14, 15]). For example of application to a human population, the turnover can
be consider as 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 occurred
during lockdown 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 livestocks
with (potential vaccinated) calf migrations between farms which influences the
epidemic.

We then complexify our analysis to analyse more finely the effect of the type
of vaccine and its efficiency on disease dynamics. The vaccine can decrease
transmission and/or favour 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 if trade-off between the vaccine efficiency to
stop transmission (infection) and disease recovery are expected. We use three
trade-off curves, and show that the convex ($\nu=\varepsilon^{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. First, we use, for illustrative purpose, Covid-19 parameters to
exemplify expected threshold for vaccination coverage for a highly
transmissible disease. Second, our model does not explicitly account for a
continuous vaccination (or a large vaccination campaigns) of individuals in a
community. Vaccination is linked in our model to the population turnover,
explaining the appearance of periodic oscillations in disease incidence (the
honey moon periods). Such periodic epidemics occur and are predicted for
Covid-19, but as a consequence of immunity waning of the various vaccines
against new variants [12]. 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.

This model contains some general conlusions which are not only applicable to
human populations, but also domesticated animals or even crops. Domesticated
animals also require vaccinations, 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 domesticated animals, the type of vaccine can be
adjusted depending on the disease, especially if it is desirable that infected
animal recover well, rather than attempting to prevent any transmission. In
addition, we also suggest that the principles of the model apply to plant
(crop) immunization. To protect plants against invasion of pathogens or pests,
one can use different biotic and synthetic chemicals to induce immunity in the
plant [29] or protect plants by spraying fungicides. 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 [30], with results mirroring our own.

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

## Authors contributions

Conception and design: HLNB, OMP, AT; Formal investigation: HLNB, OMP;
Numerical simulations: HLNB; Writing first draft: HLNB; Supervision: OMP, AT,
JMN; Revision of draft: AT, OMP, JMN.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Acknowledgements

HLB was funded by a grant from the African Institute for Mathematical
Sciences, www.nexteinstein.org, with financial support from the Government of
Canada, provided through Global Affairs Canada,
www.international.gc.ca, and the International Development Research Centre,
www.idrc.ca. OPM acknowledge the Financial support from the Alexander von
Humboldt Foundation, under the programme financed by the German Federal
Ministry of Education and Research entitled German Research Chair No
01DG15010. AT acknowledges support from the Deutsche Forschungsgemeinschaft
(DFG) through the TUM International Graduate School of Science and Engineering
(IGSSE), GSC 81, within the project GENOMIE-QADOP, and the TUM Global
Incentive Fund (exchange grant with Ghana).

## References

  * [1]
  * [2]
  * [3]
  * [4] R. Anderson, R. May, Infectious diseases of humans: dynamics and control, Oxford university press, 1991.
  * [5] P. Vale, A. Fenton, S. Brown, Limiting damage during infection: lessons from infection tolerance for novel therapeutics, PLoS biology 12 (1) (2014) e1001769.
  * [6] S. Gandon, M. Mackinnon, S. Nee, A. Read, Imperfect vaccination: some epidemiological and evolutionary consequences, Proceedings of the Royal Society of London. Series B: Biological Sciences 270 (1520) (2003) 1129–1136.
  * [7] J. Hwang, T. Zhang, A. Wang, Z. Li, Covid-19 vaccines for patients with cancer: benefits likely outweigh risks, Journal of hematology & oncology 14 (1) (2021) 1–11.
  * [8] J. Ioannidis, Benefit of covid-19 vaccination accounting for potential risk compensation, npj Vaccines 6 (1) (2021) 1–5.
  * [9] N. Dagan, N. Barda, E. Kepten, O. Miron, S. Perchik, M. Katz, M. Hernán, M. Lipsitch, B. Reis, R. Balicer, Bnt162b2 mrna covid-19 vaccine in a nationwide mass vaccination setting, New England Journal of Medicine 384 (2021) 1412–23.
  * [10] Y. G. Djatcha, S. Bowong, D. E. Houpa, J. Kurths, Mathematical analysis of the dynamical transmission of neisseria meningitidis serogroup a, International Journal of Computer Mathematics 94 (12) (2017) 2409–2434.
  * [11] B. Ashby, B. Alex, Herd immunity, Current Biology 31 (4) (2021) R174–R177.
  * [12] M. Mancuso, S. E. Eikenberry, A. B. Gumel, Will vaccine-derived protective immunity curtail covid-19 variants in the us ?, Infectious Disease Modelling 6 (2021) 1110–1134.
  * [13] A. Scherer, A. McLean, Mathematical models of vaccination, British medical bulletin 62 (1) (2002) 187–199.
  * [14] J. R. Pulliam, J. G. Dushoff, S. A. Levin, A. P. Dobson, Epidemic enhancement in partially immune populations, PLoS One 2 2 (1) (2007) e165.
  * [15] J. Knight, S. Baral, S. Schwartz, L. Wang, H. Ma, K. Young, H. Hausler, S. Mishra, Contribution of high risk groups’ unmet needs may be underestimated in epidemic models without risk turnover: A mechanistic modelling analysis, Infectious Disease Modelling 5 (2020) 549–562.
  * [16] S. Gandon, T. Day, Evidences of parasite evolution after vaccination, Vaccine 26 (2008) C4–C7.
  * [17] S. Alizon, A. Hurford, N. Mideo, M. Van Baalen, Virulence evolution and the trade‐off hypothesis: history, current state of affairs and the future, Journal of evolutionary biology 22 (2) (2009) 245–259.
  * [18] C. Cressler, D. McLeod, C. Rozins, J. Van Den Hoogen, T. Day, The adaptive evolution of virulence: a review of theoretical predictions and empirical tests, Parasitology 143 (7) (2016) 915–930.
  * [19] D. H. Anderson, Compartmental modeling and tracer kinetics, Vol. 50, Springer Science & Business Media, 2013.
  * [20] S. Gandon, T. Day, The evolutionary epidemiology of vaccination, Journal of the Royal Society Interface 4 (16) (2007) 803–817.
  * [21] J. A. Jacquez, C. P. Simon, Qualitative theory of compartmental systems, Siam Review 35 (1) (1993) 43–79.
  * [22] P. Van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Mathematical biosciences 380 (1-2) (2002) 29–48.
  * [23] Z. Shuai, P. van den Driessche, Global stability of infectious disease models using lyapunov functions, SIAM Journal on Applied Mathematics 73 (4) (2013) 1513–1532.
  * [24] H. Hethcote, H. Thieme, Stability of the endemic equilibrium in epidemic models with subpopulations, Mathematical biosciences 75 (2) (2017) 205–227.
  * [25] S. Marino, I. B. Hogue, C. J. Ray, D. E. Kirschner, A methodology for performing global uncertainty and sensitivity analysis in systems biology, Journal of theoretical biology 254 (1) (2008) 178–196.
  * [26] A. B. Gumel, C. C. McCluskey, J. Watmough, An sveir model for assessing potential impact of an imperfect anti-sars vaccine, Mathematical Biosciences and Engineering 3 (3) (2006) 485.
  * [27] A. Scherer, A. McLean, Mathematical models of vaccination, British medical bulletin 62 (1) (2002) 187–99.
  * [28] S. L. Nuismer, B. M. Althouse, R. May, J. J. Bull, S. P. Stromberg, R. Antia, Eradicating infectious disease using weakly transmissible vaccines, Proceedings of the Royal Society B 283 (1841) (2016) 20161903.
  * [29] Y. T. Dyakov, V. G. Dzhavakhiya, T. Korpela, Molecular basis of plant immunization, Elsevier, 2007.
  * [30] K. Rock, S. Brand, J. Moir, M. Keeling, Dynamics of infectious diseases, Reports on Progress in Physics 77 (2) (2014) 026602.
  * [31] J. LaSalle, Stability by Lyapunov’s direct method with applications, Vol. 50, Academic Press, 1961.
  * [32] J. P. La Salle, The stability of dynamical systems, Vol. 25, SIAM, 1976.
  * [33] H. I. Freedman, S. Ruan, M. Tang, Global stability of infectious disease models using lyapunov functions, Journal of Dynamics and Differential Equations 6 (4) (1994) 583–600.
  * [34] M. Y. Li, J. R. Graef, L. Wang, J. Karsai, Global dynamics of a seir model with varying total population size, Mathematical biosciences 160 (2) (1999) 191–213.

##  Appendix A Proof of Theorem 3.4

###### Proof.

The system (1) can be written as:

$\begin{array}[]{rl}\dfrac{\mathrm{d}a}{\mathrm{d}t}&=(F-V)a-f(a,b),\\\ \\\
\dfrac{\mathrm{d}b}{\mathrm{d}t}&=g(a,b),\end{array}$ (10)

where $a=(I_{1},I_{2})^{T}$ is the vector representing the infected classes,
$b=(S_{1},S_{2},R)^{T}$ is the vector representing the uninfected classes, the
matrices $F$ and $V$ are given as in Equation (4) and

$f(a,b)=\begin{bmatrix}\beta_{11}\Bigg{(}\dfrac{S_{1}^{0}}{N^{0}}-\dfrac{S_{1}}{N}\Bigg{)}I_{1}+\beta_{12}\Bigg{(}\dfrac{S_{1}^{0}}{N^{0}}-\dfrac{S_{1}}{N}\Bigg{)}I_{2}\\\
\\\
\beta_{21}\Bigg{(}\dfrac{S_{2}^{0}}{N^{0}}-\dfrac{S_{2}}{N}\Bigg{)}I_{1}+\beta_{22}\Bigg{(}\dfrac{S_{2}^{0}}{N^{0}}-\dfrac{S_{2}}{N}\Bigg{)}I_{2}\end{bmatrix}\hskip
5.69046pt\text{and}\hskip
5.69046ptg(a,b)=\begin{bmatrix}\theta(1-p)+\lambda_{1}S_{1}-\mu S_{1}\\\ \\\
\theta p+\lambda_{2}S_{2}-\mu S_{2}\\\ \\\ \gamma_{1}I_{1}+\gamma_{2}I_{2}-\mu
R\end{bmatrix}.$

Then,

$V^{-1}F=\begin{bmatrix}\dfrac{\beta_{11}S_{1}^{0}}{N^{0}(\mu+\gamma_{1}+d_{1})}&\dfrac{\beta_{12}S_{1}^{0}}{N^{0}(\mu+\gamma_{1}+d_{1})}\\\
\\\
\dfrac{\beta_{21}S_{2}^{0}}{N^{0}(\mu+\gamma_{2}+d_{2})}&\dfrac{\beta_{22}S_{2}^{0}}{N^{0}(\mu+\gamma_{2}+d_{2})}\end{bmatrix},$

and the left eigenvector of $V^{-1}F$, $(\omega_{1},\omega_{2})$ associated
with the eigenvalue $\mathcal{R}_{0}$ is given by:

$\omega_{1}=1$ and
$\omega_{2}=\dfrac{N^{0}(\mu+\gamma_{2}+d_{2})}{\beta_{21}S_{2}^{0}}\Bigg{(}\mathcal{R}_{0}-\dfrac{\beta_{11}S_{1}^{0}}{N^{0}(\mu+\gamma_{1}+d_{1})}\Bigg{)}$
since

$(\omega_{1},\omega_{2})V^{-1}F=\mathcal{R}_{0}(\omega_{1},\omega_{2}).$

Let us consider the following Lyapunov function:

$\displaystyle Q=$
$\displaystyle(\omega_{1},\omega_{2})V^{-1}(I_{1},I_{2})^{T}$ $\displaystyle=$
$\displaystyle\dfrac{I_{1}}{\mu+\gamma_{1}+d_{1}}+\Bigg{(}\mathcal{R}_{0}-\dfrac{\beta_{11}S_{1}^{0}}{N^{0}(\mu+\gamma_{1}+d_{1})}\Bigg{)}\dfrac{N^{0}I_{2}}{\beta_{21}S_{2}^{0}}.$
(11)

Then the derivative of $Q$ with respect to $t$ yields,

$Q^{{}^{\prime}}=(\mathcal{R}_{0}-1)(\omega_{1},\omega_{2})^{T}a-(\omega_{1},\omega_{2})^{T}V^{-1}f(a,b).$

Since $(\omega_{1},\omega_{2})\geqslant 0,$ $V^{-1}\geqslant 0$ and
$f(a,b)\geqslant 0$ in $\Omega$, then
$(\omega_{1},\omega_{2})^{T}V^{-1}f(a,b)\geqslant 0.$ Therefore,
$Q^{{}^{\prime}}\leqslant 0$ in $\Omega$ if $\mathcal{R}_{0}\leqslant 1$ and
$Q$ is a Lyapunov function for the system (1). By LaSalle’s invariance
principle [31, 32], $Q^{0}$ is GAS in $\Omega$.

If $\mathcal{R}_{0}>1,$ then
$Q^{{}^{\prime}}=(\mathcal{R}_{0}-1)(\omega_{1},\omega_{2})^{T}a>0$ provided
that $a>0$ and $b=(S_{1}^{0},S_{2}^{0},0)$. By continuity, $Q^{{}^{\prime}}>0$
in the neighborhood of $Q^{0}$. Solutions in positive cone sufficiently close
to $Q^{0}$ move away from $Q^{0}$, implying that $Q^{0}$ is unstable. Thus,
the model system (1) is uniformly persistent [33, 34]. Uniform persistence and
the positively invariance of $\Omega$ imply the existence of an endemic
equilibrium. ∎

##  Appendix B Proof of Theorem 3.5

###### 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\Bigg{(}\dfrac{S_{1}}{S_{1}^{*}}\Bigg{)}$,
$L_{2}=S_{2}-S_{2}^{*}-S_{2}^{*}\log\Bigg{(}\dfrac{S_{2}}{S_{2}^{*}}\Bigg{)}$,
$L_{3}=I_{1}-I_{1}^{*}-I_{1}^{*}\log\Bigg{(}\dfrac{S_{3}}{S_{3}^{*}}\Bigg{)}$
and
$L_{4}=I_{2}-I_{2}^{*}-I_{2}^{*}\log\Bigg{(}\dfrac{I_{4}}{I_{4}^{*}}\Bigg{)}$.

Using the inequality $1-z+\log(z)\leqslant 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}[]{llll}L_{1}^{{}^{\prime}}&=\Bigg{(}1-\dfrac{S^{*}_{1}}{S_{1}}\Bigg{)}\dfrac{\mathrm{d}S_{1}}{\mathrm{d}t}\\\
\\\
&=\Bigg{(}1-\dfrac{S^{*}_{1}}{S_{1}}\Bigg{)}\Bigg{[}\beta_{11}\dfrac{S_{1}^{*}I_{1}^{*}}{N^{*}}-\beta_{11}\dfrac{S_{1}I_{1}}{N}+\beta_{12}\dfrac{S_{1}^{*}I_{2}^{*}}{N^{*}}-\beta_{12}\dfrac{S_{1}I_{2}}{N}-\mu
S_{1}+\mu S^{*}_{1}\Bigg{]}\\\ \\\
&=-\dfrac{\mu(S_{1}-S^{*}_{1})^{2}}{S_{1}}+\beta_{11}\dfrac{S_{1}^{*}I_{1}^{*}}{N^{*}}\Bigg{[}1-\dfrac{S^{*}_{1}}{S_{1}}-\dfrac{S_{1}I_{1}N^{*}}{S_{1}^{*}I_{1}^{*}N}+\dfrac{I_{1}N^{*}}{I_{1}^{*}N}\Bigg{]}+\beta_{12}\dfrac{S_{1}^{*}I_{2}^{*}}{N^{*}}\Bigg{[}1-\dfrac{S^{*}_{1}}{S_{1}}-\dfrac{S_{1}I_{2}N^{*}}{S_{1}^{*}I_{2}^{*}N}+\dfrac{I_{2}N^{*}}{I_{2}^{*}N}\Bigg{]}.\\\
\\\ \text{Then}\quad
L_{1}^{{}^{\prime}}\leqslant&\beta_{11}\dfrac{S_{1}^{*}I_{1}^{*}}{N^{*}}\Bigg{[}\dfrac{I_{1}N^{*}}{I_{1}^{*}N}-\log\Bigg{(}\dfrac{I_{1}N^{*}}{I_{1}^{*}N}\Bigg{)}-\dfrac{S_{1}I_{1}N^{*}}{S_{1}^{*}I_{1}^{*}N}+\log\Bigg{(}\dfrac{S_{1}I_{1}N^{*}}{S_{1}^{*}I_{1}^{*}N}\Bigg{)}\Bigg{]}\\\
\\\
&+\beta_{12}\dfrac{S_{1}^{*}I_{2}^{*}}{N^{*}}\Bigg{[}\dfrac{I_{2}N^{*}}{I_{2}^{*}N}-\log\Bigg{(}\dfrac{I_{2}N^{*}}{I_{2}^{*}N}\Bigg{)}-\dfrac{S_{1}I_{2}N^{*}}{S_{1}^{*}I_{2}^{*}N}+\log\Bigg{(}\dfrac{S_{1}I_{2}N^{*}}{S_{1}^{*}I_{2}^{*}N}\Bigg{)}\Bigg{]}.\end{array}$
(12)

We can also deduce in an analogous way:

$\begin{array}[]{llll}L_{2}^{{}^{\prime}}&\leqslant\beta_{22}\dfrac{S_{2}^{*}I_{2}^{*}}{N^{*}}\Bigg{[}\dfrac{I_{2}N^{*}}{I_{2}^{*}N}-\log\Bigg{(}\dfrac{I_{2}N^{*}}{I_{2}^{*}N}\Bigg{)}-\dfrac{S_{2}I_{2}N^{*}}{S_{2}^{*}I_{2}^{*}N}+\log\Bigg{(}\dfrac{S_{2}I_{2}N^{*}}{S_{2}^{*}I_{2}^{*}N}\Bigg{)}\Bigg{]}\\\
\\\
&+\beta_{21}\dfrac{S_{2}^{*}I_{1}^{*}}{N^{*}}\Bigg{[}\dfrac{I_{1}N^{*}}{I_{1}^{*}N}-\log\Bigg{(}\dfrac{I_{1}N^{*}}{I_{1}^{*}N}\Bigg{)}-\dfrac{S_{2}I_{1}N^{*}}{S_{2}^{*}I_{1}^{*}N}+\log\Bigg{(}\dfrac{S_{2}I_{1}N^{*}}{S_{2}^{*}I_{1}^{*}N}\Bigg{)}\Bigg{]}.\end{array}$
(13)

We also have

$\begin{array}[]{llll}L_{3}^{{}^{\prime}}&=\Bigg{(}1-\dfrac{I^{*}_{1}}{I_{1}}\Bigg{)}\dfrac{\mathrm{d}I_{1}}{\mathrm{d}t}\\\
\\\
&=\Bigg{(}1-\dfrac{I^{*}_{1}}{I_{1}}\Bigg{)}\Bigg{[}\beta_{11}\dfrac{S_{1}I_{1}}{N}+\beta_{12}\dfrac{S_{1}I_{2}}{N}-(\mu+\gamma_{1}+d_{1})I_{1}\Bigg{]}\\\
\\\
&=\Bigg{(}1-\dfrac{I^{*}_{1}}{I_{1}}\Bigg{)}\Bigg{[}\beta_{11}\dfrac{S_{1}I_{1}}{N}+\beta_{12}\dfrac{S_{1}I_{2}}{N}-\beta_{11}\dfrac{S_{1}^{*}I_{1}}{N^{*}}+\beta_{12}\dfrac{S_{1}^{*}I_{2}^{*}I_{1}}{N^{*}I_{1}^{*}}\Bigg{]}\\\
\\\
&=\beta_{11}\dfrac{S_{1}^{*}I_{1}^{*}}{N^{*}}\Bigg{[}\dfrac{S_{1}I_{1}N^{*}}{S_{1}^{*}I_{1}^{*}N}-\dfrac{S_{1}N^{*}}{S_{1}^{*}N}-\dfrac{I_{1}}{I_{1}^{*}}+1\Bigg{]}+\beta_{12}\dfrac{S_{1}^{*}I_{2}^{*}}{N^{*}}\Bigg{[}\dfrac{S_{1}I_{2}N^{*}}{S_{1}^{*}I_{2}^{*}N}-\dfrac{S_{1}I_{1}^{*}I_{2}N^{*}}{S_{1}^{*}I_{1}I_{2}^{*}N}-\dfrac{I_{1}}{I_{1}^{*}}+1\Bigg{]},\\\
\\\
L_{3}^{{}^{\prime}}&\leqslant\beta_{11}\dfrac{S_{1}^{*}I_{1}^{*}}{N^{*}}\Bigg{[}\dfrac{S_{1}I_{1}N^{*}}{S_{1}^{*}I_{1}^{*}N}-\log\Bigg{(}\dfrac{S_{1}I_{1}N^{*}}{S_{1}^{*}I_{1}^{*}N}\Bigg{)}-\dfrac{I_{1}}{I_{1}^{*}}+\log\Bigg{(}\dfrac{I_{1}}{I_{1}^{*}}\Bigg{)}\Bigg{]}\\\
\\\
&+\beta_{12}\dfrac{S_{1}^{*}I_{2}^{*}}{N^{*}}\Bigg{[}\dfrac{S_{1}I_{2}N^{*}}{S_{1}^{*}I_{2}^{*}N}-\log\Bigg{(}\dfrac{S_{1}I_{2}N^{*}}{S_{1}^{*}I_{2}^{*}N}\Bigg{)}-\dfrac{I_{1}}{I_{1}^{*}}+\log\Bigg{(}\dfrac{I_{1}}{I_{1}^{*}}\Bigg{)}\Bigg{]}.\end{array}$
(14)

Similarly, we obtain

$\begin{array}[]{llll}L_{4}^{{}^{\prime}}&\leqslant\beta_{22}\dfrac{S_{2}^{*}I_{2}^{*}}{N^{*}}\Bigg{[}\dfrac{S_{2}I_{2}N^{*}}{S_{2}^{*}I_{2}^{*}N}-\log\Bigg{(}\dfrac{S_{2}I_{2}N^{*}}{S_{2}^{*}I_{2}^{*}N}\Bigg{)}-\dfrac{I_{2}}{I_{2}^{*}}+ln\dfrac{I_{2}}{I_{2}^{*}}\Bigg{]}\\\
\\\
&+\beta_{21}\dfrac{S_{2}^{*}I_{1}^{*}}{N^{*}}\Bigg{[}\dfrac{S_{2}I_{1}N^{*}}{S_{2}^{*}I_{1}^{*}N}-\log\Bigg{(}\dfrac{S_{2}I_{1}N^{*}}{S_{2}^{*}I_{1}^{*}N}\Bigg{)}-\dfrac{I_{2}}{I_{2}^{*}}+\log\Bigg{(}\dfrac{I_{2}}{I_{2}^{*}}\Bigg{)}\Bigg{]}.\end{array}$
(15)

Therefore, by adding (12), (13), (14) and (15) we deduce

$\displaystyle L^{{}^{\prime}}\leqslant$
$\displaystyle\Bigg{(}-\dfrac{I_{1}N^{*}}{I_{1}^{*}N}+\log\Bigg{(}\dfrac{I_{1}N^{*}}{I_{1}^{*}N}\Bigg{)}\Bigg{)}\Bigg{(}-\beta_{11}\dfrac{S_{1}^{*}I_{1}^{*}}{N^{*}}-\beta_{21}\dfrac{S_{2}^{*}I_{1}^{*}}{N^{*}}\Bigg{)}$
$\displaystyle+\Bigg{(}-\dfrac{I_{2}N^{*}}{I_{2}^{*}N}+\log\Bigg{(}\dfrac{I_{2}N^{*}}{I_{2}^{*}N}\Bigg{)}\Bigg{)}\Bigg{(}-\beta_{12}\dfrac{S_{1}^{*}I_{2}^{*}}{N^{*}}-\beta_{22}\dfrac{S_{2}^{*}I_{2}^{*}}{N^{*}}\Bigg{)}$
$\displaystyle+\Bigg{(}-\dfrac{I_{1}}{I_{1}^{*}}+\log\Bigg{(}\dfrac{I_{1}}{I_{1}^{*}}\Bigg{)}\Bigg{)}\Bigg{(}\beta_{11}\dfrac{S_{1}^{*}I_{1}^{*}}{N^{*}}+\beta_{12}\dfrac{S_{1}^{*}I_{2}^{*}}{N^{*}}\Bigg{)}$
$\displaystyle+\Bigg{(}-\dfrac{I_{2}}{I_{2}^{*}}+\log\Bigg{(}\dfrac{I_{2}}{I_{2}^{*}}\Bigg{)}\Bigg{)}\Bigg{(}\beta_{22}\dfrac{S_{2}^{*}I_{2}^{*}}{N^{*}}+\beta_{21}\dfrac{S_{2}^{*}I_{1}^{*}}{N^{*}}\Bigg{)}.$
$\displaystyle\text{Then}\quad L^{{}^{\prime}}\leqslant$ $\displaystyle
0,\hskip 5.69046pt\text{since}\hskip 5.69046pt-z+\log(z)\leqslant-1,\hskip
5.69046pt\forall z>0.$

Since $\\{Q^{*}\\}$ is the only invariant subset in $\Omega$ where $L=0$,
therefore by LaSalle’s invariance principle [32], $Q^{*}$ is GAS in $\Omega$.
∎

##  Appendix C Tables

Table 3: 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: $I_{1}+I_{2}$
---|---|---
| Min | Baseline | Max |  $t=50$ days |  $t=100$ days |  $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 4: 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: $I_{1}+I_{2}$
---|---|---
| Min | Baseline | Max |  $t=50$ days |  $t=100$ days |  $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 5: 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: $I_{1}+I_{2}$
---|---|---
| Min | Baseline | Max |  $t=50$ days |  $t=100$ days |  $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 6: 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: $I_{1}+I_{2}$
---|---|---
| Min | Baseline | Max |  $t=50$ days |  $t=100$ days |  $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$

##  Appendix D Figures

Figure 7: Scatter plots of $\mathcal{R}_{0}$ with a strong turnover as a
function of $\varepsilon$, $\nu$ and $p$ .The parameters are $\theta=1000$,
$\beta_{11}=0.35$, $\beta_{12}=0.28$, $\beta_{21}=0.175,\beta_{22}=0.14$,
$\mu=0.09$, $d_{1}=0.0008$, $d_{2}=0.0001$, $\gamma_{1}=0.065$,
$\gamma_{2}=0.13$. Figure 8: Slice planes of $\mathcal{R}_{0}$ orthogonal to
the p-axis at the values 0.2, 0.5, 0.8 with a strong turnover. The parameters
are $\theta=1000$, $\beta_{11}=0.35$, $\beta_{12}=0.28$,
$\beta_{21}=0.175,\beta_{22}=0.14$, $\mu=0.09$, $d_{1}=0.0008$,
$d_{2}=0.0001$, $\gamma_{1}=0.065$, $\gamma_{2}=0.13$.

(a)

(b)

(c)

Figure 9: Contour plots of the basic reproduction number ($\mathcal{R}_{0}$)
of the model (1) with a weak turnover as a function of: (a) vaccine coverage,
$p$, and vaccine efficiency on the transmission, $\varepsilon$ (fixed
$\nu=0.5$); (b) vaccine coverage, $p$, and vaccine efficiency on the ability
to enhance recovery, $\nu$ (fixed $\varepsilon=0.5$); (c) vaccine efficiency
on the ability of being recovered, $\nu$, and vaccine efficiency on the
transmission, $\varepsilon$ (fixed $p=0.5$). The parameters are $\theta=10$,
$\beta_{11}=0.35$, $\beta_{12}=0.28$, $\beta_{21}=0.175,\beta_{22}=0.14$,
$\mu=0.0009$, $d_{1}=0.0008$, $d_{2}=0.0001$, $\gamma_{1}=0.065$,
$\gamma_{2}=0.13$ Figure 10: Slice planes of $\mathcal{R}_{0}$ orthogonal to
the p-axis at the values 0.2, 0.5, 0.8 with a weak turnover. The parameters
are $\theta=10$, $\beta_{11}=0.35$, $\beta_{12}=0.28$,
$\beta_{21}=0.175,\beta_{22}=0.14$, $\mu=0.0009$, $d_{1}=0.0008$,
$d_{2}=0.0001$, $\gamma_{1}=0.065$, $\gamma_{2}=0.13$. Figure 11: Scatter
plots of $\mathcal{R}_{0}$ with a weak turnover as a function of
$\varepsilon$, $\nu$ and $p$ .The parameters are $\theta=10$,
$\beta_{11}=0.35$, $\beta_{12}=0.28$, $\beta_{21}=0.175,\beta_{22}=0.14$,
$\mu=0.0009$, $d_{1}=0.0008$, $d_{2}=0.0001$, $\gamma_{1}=0.065$,
$\gamma_{2}=0.13$.

(a)

(b)

(c)

Figure 12: Contour plot of the basic reproduction number ($\mathcal{R}_{0}$)
of the model (1) with a weak turnover as a function of vaccine coverage, $p$,
and vaccine efficiency on the transmission, $\varepsilon$ when: (a)
$\nu=\varepsilon^{2}$(convex relationship); (b)
$\nu=\sqrt{\varepsilon}$(concave relationship); (c)$\nu=\varepsilon$(linear
relationship). The parameters are $\theta=1000$, $\beta_{11}=0.35$,
$\beta_{12}=0.28$, $\beta_{21}=0.175,\beta_{22}=0.14$, $\mu=0.09$, $p=0.5$,
$d_{1}=0.0008$, $d_{2}=0.0001$, $\gamma_{1}=0.065$, $\gamma_{2}=0.13$.