††footnotetext: MSC2020: 00A71, 34D20, 37M05, 37N25, 92D30.

# Dynamics of a mathematical model of virus spreading incorporating the effect
of a vaccine

Aytül Gökçe1 and Burcu Gürbüz∗,2 and Alan D. Rendall3 A. Gökçe1 Ordu
University, Faculty of Science and Letters, Department of Mathematics, 52200,
Ordu, Turkey<EMAIL_ADDRESS>B. Gürbüz∗,2 Institut für Mathematik,
Johannes Gutenberg-Universität, Staudingerweg 9, 55099, Mainz, Germany
<EMAIL_ADDRESS>A. D. Rendall3 Institut für Mathematik, Johannes
Gutenberg-Universität, Staudingerweg 9, 55099, Mainz, Germany  rendall@uni-
mainz.de

(Date: September 3, 2024)

###### Abstract.

The COVID-19 pandemic led to widespread interest in epidemiological models. In
this context the role of vaccination in influencing the spreading of the
disease is of particular interest. There has also been a lot of debate on the
role of non-pharmaceutical interventions such as the disinfection of surfaces.
We investigate a mathematical model for the spread of a disease which includes
both imperfect vaccination and infection due to virus in the environment. The
latter is studied with the help of two phenomenological models for the force
of infection. In one of these models we find that backward bifurcations take
place so that for some parameter values an endemic steady state exists
although the basic reproduction ratio $R_{0}$ is less than one. We also prove
that in that case there can exist more than one endemic steady state. In the
other model all generic transcritical bifurcations are forward bifurcations so
that these effects cannot occur. Thus we see that the occurrence of backward
bifurcations, which can be important for disease control strategies, is
dependent on the details of the function describing the force of infection. By
means of simulations the predictions of this model are compared with data for
COVID-19 from Turkey. A sensitivity analysis is also carried out.

1\. Introduction

Since the beginning of epidemiology mathematical models have played a central
role. This can be seen in the groundbreaking work of Ronald Ross and Hilda
Hudson on the eradication of malaria [24], [25], [26]. In that work the
authors identified a threshold for the persistence of a disease which can be
seen as the ancestor of the basic reproduction ratio $R_{0}$ which is so
important in epidemiology today. The COVID-19 epidemic caused a surge of work
where epidemiological models were defined, simulated and subjected to rigorous
mathematical analysis. Due to the urgency of the situation this development
took place in a rather disorganized way. Now it is time to consolidate and
extend the things learned at that time so as to be prepared as well as
possible for future epidemics.

In this paper we study a model for the spread of an infectious disease in a
human population which includes an imperfect vaccination and takes into
account infections due to virus particles in the environment. In particular we
are thinking of fomites, objects in the environment which are contaminated
with virus and which are not humans or animals. Here we may think of the
contamination by hands touching doorknobs [29] or infections spreading in
hospitals [27]. The question of how important this route of infection is for
COVID-19 has been a subject of much discussion. The consensus appears to be
that it is of secondary importance but this may be different for other
diseases [2, 20].

In Section 2 we define the model which is of central interest in this paper
and establish some basic properties of its solutions. The model contains a
response function which describes how the concentration of virus in the
environment affects the rate of infection by this route. This function depends
on an integer $n\geq 1$. The motivation for the choice of this function is
also discussed. In Section 2.1 it is shown that this model has a unique
disease-free steady state (Lemma 3.1) and the stability of that state is
determined using the next generation matrix. The basic reproduction ratio
$R_{0}$ is computed for this model and it is shown to be a decreasing function
of the vaccination rate. The existence of backward bifurcations is analysed
using the method of van den Driessche and Watmough [12]. It is proved that
generic backward bifurcations occur in the case $n=2$ but not otherwise
(Theorem 3.1). In particular they do not occur in the case $n=1$, where the
function describing the force of infection is one which had previously been
considered in the literature [6]. Note that an existing model for imperfect
vaccination [15] also exhibits no backward bifurcations. A backward
bifurcation is often accompanied by the presence of more than one positive
steady state for given values of the parameters (cf. [19]). In many cases of
backward bifurcations simulations show not only that for a certain choice of
parameters positive steady states exist for $R_{0}<1$ but also that two
positive steady states can occur. It is proved in Section 4 that there are
parameters for which our model with $n=2$ exhibits the latter behaviour. In
Section 5 it is shown that solutions of the model can be fitted to COVID-19
data from Turkey. Section 6 carries out a sensitivity analysis of the model.

2\. The model

The model considered in what follows is a generalization of one introduced in
[15] to study the effects of vaccination against SARS and is given by the
following equations:

$\displaystyle\frac{dS}{dt}$ $\displaystyle=\Lambda-\beta SI-\sigma
S+(1-\lambda)t^{\prime}V-\alpha_{1}Sg(C,\kappa)-\mu S,$ (1)
$\displaystyle\frac{dE}{dt}$ $\displaystyle=\beta SI+\epsilon IV-\xi E-\mu
E+\alpha_{1}Sg(C,\kappa)+\alpha_{2}Vg(C,\kappa),$ (2)
$\displaystyle\frac{dI}{dt}$ $\displaystyle=\xi E-\delta I-dI-\mu I,$ (3)
$\displaystyle\frac{dV}{dt}$ $\displaystyle=\sigma S-\epsilon
IV-(t^{\prime}+\mu)V-\alpha_{2}Vg(C,\kappa),$ (4) $\displaystyle\frac{dR}{dt}$
$\displaystyle=\delta I-\mu R+\lambda t^{\prime}V,$ (5)
$\displaystyle\frac{dC}{dt}$ $\displaystyle=\varphi I-\omega C,$ (6)

where

$g(C,\kappa)=\frac{C^{n}}{C^{n}+\kappa}.$

The meaning of the parameters in this model is described in Table 1. The model
of [15] was called an SVEIR model after the names of its five unknowns $S$,
$V$, $E$, $I$ and $R$. These are the numbers of susceptible, vaccinated,
exposed, infectious and recovered individuals, respectively. We augment this
by an additional variable $C$ representing the concentration of the virus in
the environment. In both models the vaccination is imperfect but the
imperfection is of a different kind. Correspondingly the class $V$ has a
different interpretation in the two models. In [15] the class $V$ consists of
individuals who have been vaccinated at some time. The effect of the
vaccination is to lower the rate at which they get infected compared to
unvaccinated individuals. In the model (1)-(6) the class $V$ consists of
individuals who have received a vaccination but where the vaccination has not
yet had time to become fully effective. After that time either the vaccination
provides complete protection or it has not been effective and the individual
returns to the susceptible class. For biological reasons the inequalities
$\epsilon\leq\beta$ and $\alpha_{2}\leq\alpha_{1}$ are assumed, which means
that vaccinated individuals are no more likely to be infected than
unvaccinated individuals, either by infected individuals or by contact with
their surroundings.

Mathematically the model of [15], up to a different notation, can be obtained
from our model by setting $t^{\prime}=\alpha_{1}=\alpha_{2}=0$ and discarding
the equation for $S$. This is possible since when the parameters just listed
are zero the equations for the first five variables to not depend on $S$. In
the model (1)-(6) the imperfection of the vaccination is expressed as follows.
Individuals leave the vaccinated state at rate $t^{\prime}$, the vaccination
having been successful with probability $\lambda$. This way of modelling an
imperfect vaccination was previously used in [3]. The other additional effect
taken account of in (1)-(6) is related to infection by virus in the
environment. It is expressed by the terms containing the factors $\alpha_{1}$
and $\alpha_{2}$ relating to the unvaccinated and vaccinated individuals,
respectively. This type of effect was included in a model of [6]. In that
paper the function $g$ written above with $n=1$ was used as a phenomenological
description of the rate of infection in this process.

It is worth taking some time to discuss the status of this type of
phenomenological description. It is used in defining response functions in
various parts of biology. In biochemistry the case $n=1$ is called a
Michaelis-Menten function while the case $n\geq 2$ is called a Hill function.
In predator-prey models in ecology the case $n=1$ is called Holling type II
while the case $n\geq 2$ is called Holling type III. Holling type I denotes a
linear response function, usually with a cut-off. A general discussion of
response functions in epidemiology is given in Chapter 10 of the book of
Diekmann and Heesterbeek [10], whereby the authors make clear from the
beginning that they do not claim to give a definitive answer to the questions
they are raising.

Suppose that there is a source of infection with intensity $Z$ and a
population $S$ of susceptibles. Let $F(Z,S)$ be the rate of infection. In
principle this could be any function. Let us suppose that $F$ depends linearly
on $S$ but initially allow its dependence on $Z$ to be arbitrary. Thus
$F(Z,S)=sf(Z)$ for some function $f$. What properties should the function $f$
have? It should be positive for $Z$ positive and zero for $Z=0$. It should be
non-decreasing. It is reasonable to assume that it is bounded. The simplest
type of function satisfying these requirements is one of the form
$f(Z)=\frac{aZ}{Z+b}$. Another situation in which a response function is of
relevance is the predation rate in a predator-prey model. There the analogue
of $I$ is the density of predators while the analogue of $Z$ is the density of
prey. In that case a function of the form just considered is called Holling
type II. In that context Holling type I is a function which is linear up to a
threshold value and then constant. Holling type III corresponds to $Z$ being
replaced by $Z^{p}$. This argument for introducing a function $f$ of this form
is purely phenomenological. Holling had a mechanistic argument to motivate his
type II function. There have also been attempts to motivate the type III
function by mechanistic considerations (cf. [9], [16]). We are not aware that
this has been done in epidemiology. The function corresponding to Holling type
II was introduced to epidemiological models by Dietz [11], without a
mechanistic background. Holling’s mechanistic approach does not apply to
epidemiological models. Diekmann and Heesterbeek [10] discuss mechanistic
approaches to the Holling type II function in epidemiology. In fact in a model
case they derive something which is not a rational function. We have not found
a paper where the Holling type III function is used in epidemiology.

Table 1.  Parameters used in the model.  parameters | biological meaning
---|---
$\Lambda$ | recruitment rate
$\beta$ | effective contact rate with $\beta S$ new susceptible individuals per unit time
$\alpha_{1}$ | transmission ratio of the virus from the environment to susceptible individuals
| that enter the exposed class
$\alpha_{2}$ | transmission ratio of the virus from the environment to vaccinated individuals
| (not fully immunised) that may enter the exposed class
$\epsilon$ | the rate at which a vaccinated individual (not fully immunised) becomes exposed
| after being in contact with an infected individual
$\delta$ | the rate of recovered individuals
$\mu$ | natural mortality rate
$\xi$ | rate of development of clinical symptoms
$d$ | disease induced fatality rate
$\sigma$ | vaccination rate of susceptible individuals (the first shot)
$\varphi$ | the virus exposure rate
$\lambda$ | the efficiency of the vaccine
$t^{\prime-1}$ | the mean amount of time that is spent in the vaccinated class before developing
| an immune response and moving to the recovery class.
$\omega$ | the rate of decay in the virus density

The right hand sides of equations (1)-(6) are smooth and hence for any initial
values at a given time they have a unique solution on some time interval.
Because of the interpretation of the unknowns we are interested in solutions
which are non-negative at all times. This is true provided the initial data
are non-negative (cf. [22], Lemma 1). Let $N(t)=S(t)+E(t)+I(t)+V(t)+R(t)$.
Then

$\frac{{\rm d}N}{{\rm d}t}=\Lambda-\mu N-dI\leq\Lambda-\mu N.$ (7)

This implies that $N$ remains bounded on any finite time interval. Hence on
such an interval all variables other than $C$ are bounded. It then follows
from (6) that $C$ is also bounded. As a consequence the solutions exist
globally in the future. Using the differential inequality for $N$ again shows
that

$\limsup\limits_{t\to\infty}N(t)\leq\frac{\Lambda}{\mu}.$ (8)

It then follows from (6) that
$\limsup_{t\to\infty}C(t)\leq\frac{\Lambda\varphi}{\mu\omega}$.

2.1. The disease-free steady state

Consider a boundary steady state $(S^{*},E^{*},I^{*},V^{*},R^{*},C^{*})$ of
the system (1)-(6), i.e. a time-independent solution for which at least one of
the unknowns is zero.

Lemma 2.1 The model (1)-(6) has a unique boundary steady state.

Proof Let $(S^{*},E^{*},I^{*},V^{*},R^{*},C^{*})$ be a boundary steady state.
If $S^{*}=0$ then (1) gives a contradiction and so $S^{*}\neq 0$. If $V^{*}=0$
then (4) implies that $S^{*}=0$. Hence in fact $V^{*}\neq 0$. If $R^{*}=0$
then (6) implies that $V^{*}=0$. Hence in fact $R^{*}\neq 0$. It follows from
the other three equations that for a steady state the equations $E^{*}=0$,
$I^{*}=0$ and $C^{*}=0$ are all equivalent to each other. Hence any boundary
steady state is of the form $(S^{*},0,0,V^{*},R^{*},0)$. In this case the
steady state equations are equivalent to the system

$\displaystyle\Lambda-\sigma S^{*}+(1-\lambda)t^{\prime}V^{*}-\mu S^{*}=0,$
$\displaystyle\sigma S^{*}-(t^{\prime}+\mu)V^{*}=0,$ (9) $\displaystyle\mu
R^{*}+\lambda t^{\prime}V^{*}=0,$

and these can be solved to give

$\displaystyle S^{*}$
$\displaystyle=\frac{\Lambda(t^{\prime}+\mu)}{\mu(\sigma+t^{\prime}+\mu)+\lambda
t^{\prime}\sigma},$ (10) $\displaystyle V^{*}$
$\displaystyle=\frac{\Lambda\sigma}{\mu(\sigma+t^{\prime}+\mu)+\lambda
t^{\prime}\sigma},$ (11) $\displaystyle R^{*}$ $\displaystyle=\frac{\lambda
t^{\prime}\Lambda\sigma}{\mu(\mu(\sigma+t^{\prime}+\mu)+\lambda
t^{\prime}\sigma)}.$ (12)

Thus there exists a unique boundary steady state of this model. $\blacksquare$

All the variables corresponding to the presence of infection are zero in this
state and so we call it the disease-free steady state.

To illustrate how a solution of this model corresponding to an epidemic might
look we show the results of a simulation with biologically motivated
parameters. Table 2 lists the references which were used either as direct
sources or guidelines for the choice of the parameters. Figure 1 demonstrates
the dynamics of populations during an epidemic for $300$ days. Susceptible
($S$), Infected ($I$) and Recovered ($R$) populations respectively associated
with blue, red and green lines. The initial conditions are chosen as
$S_{0}=61098000,V_{0}=18500000,E_{0}=2200000,I_{0}=1200000$ and $R_{0}=2000$
and parameters are given in Table 2. It should be noted that these initial
conditions and parameters are only selected for illustrative purposes and may
not be epidemiologically realistic.

Table 2. Biologically meaningful parameters used in Fig. 1. Parameters | Value ($n=1$) | Unit | Source
---|---|---|---
$\Lambda$ | $3032$ | day-1 | assumed based on [5, 14, 15]
$\beta$ | $0.15\times 10^{-8}$ | day-1 | assumed based on [15]
$\mu$ | $3.653\times 10^{-5}$ | day-1 | assumed based on [6, 15]
$\epsilon$ | $0.15\times 10^{-8}$ | day-1 | assumed based on [3, 15]
$t^{\prime}$ | $1/120$ | day-1 | assumed based on [3]
$\lambda$ | $0.8$ | day-1 | assumed based on [3, 4, 21]
$d$ | $0.02$ | day-1 | assumed based on [6, 23]
$\alpha_{1}$ | $0.01$ | day-1 | assumed
$\alpha_{2}$ | $0.01$ | day-1 | assumed
$\omega$ | $4$ | day-1 | [6]
$\varphi$ | $2$ | day-1 | [6]
$\kappa$ | $20000$ | copies $/$day | assumed based on [6]
$\xi$ | $0.125$ | day-1 | [15]
$\delta$ | $0.06$ | day-1 | assumed based on [3, 15]
$\sigma$ | $0.01$ | day-1 | [3]

Here a dramatic increase can be seen in the number of infected individuals
until day $45$, then a gentle decline appears for the population of infected
individuals. As is observed from the graph the number of susceptible
individuals slowly decreases to $16000000$ and while the number of recovered
individuals rises above $62000000$ at day $290$. The total population is taken
as $83$ million.

Figure 1. Simulation result of the model (1)-(6) with initial data and
parameters given in the text.

3\. Stability of the disease-free steady state

Linearisation around the disease-free steady state leads to the Jacobian
matrix

$\mathcal{J}=[J_{ij}]_{6\times 6}\biggr{\rvert}_{E^{*}},\ \mbox{for}\
i,j=1,2,...,6,$ (13)

where $E^{*}=(S^{*},0,0,V^{*},R^{*},0)$ and

$\displaystyle J_{11}$ $\displaystyle=-\sigma-\mu,\ J_{12}=0,\ J_{13}=-\beta
S,\ J_{14}=(1-\lambda)t^{\prime},\ J_{15}=0,\
J_{16}=-{\alpha_{1}S}\delta_{1n}/{\kappa},$ $\displaystyle J_{21}$
$\displaystyle=0,\ J_{22}=-\xi-\mu,\ J_{23}=\epsilon V+\beta S,\ J_{24}=0,\
J_{25}=0,\ J_{26}=(\alpha_{1}S+\alpha_{2}V)\delta_{1n}/\kappa,$ $\displaystyle
J_{31}$ $\displaystyle=0,\ J_{32}=\xi,\ J_{33}=-(\delta+d+\mu),\ J_{34}=0,\
J_{35}=0,\ J_{36}=0,$ $\displaystyle J_{41}$ $\displaystyle=\sigma,\
J_{42}=0,\ J_{43}=-\epsilon V,\ J_{44}=-(t^{\prime}+\mu),\ J_{45}=0,\
J_{46}=-\alpha_{2}V\delta_{1n}/\kappa,$ $\displaystyle J_{51}$
$\displaystyle=0,\ J_{52}=0,\ J_{53}=\delta,\ J_{54}=\lambda t^{\prime},\
J_{55}=-\mu,J_{56}=0,$ $\displaystyle J_{61}$ $\displaystyle=0,\ J_{62}=0,\
J_{63}=\varphi,\ J_{64}=0,\ J_{65}=0,\ J_{66}=-\omega.$

Here $\delta_{1n}$ is a Kronecker delta. Then we have the matrix

$\displaystyle\mathcal{J}$
$\displaystyle=\begin{bmatrix}J_{11}&0&J_{13}&J_{14}&0&J_{16}\\\
0&J_{22}&J_{23}&0&0&J_{26}\\\ 0&J_{32}&J_{33}&0&0&0\\\
J_{41}&0&J_{43}&J_{44}&0&J_{46}\\\ 0&0&J_{53}&J_{54}&J_{55}&0\\\
0&0&J_{63}&0&0&J_{66}\\\ \end{bmatrix}$

It is clear that one of the eigenvalues of the Jacobian is $J_{55}$. Moreover,
removing the fifth row and column and interchanging the second row and column
with the fourth leads to a matrix with block diagonal structure. Thus two
further eigenvalues of the Jacobian can be obtained as the eigenvalues of the
matrix

$\displaystyle\overline{\mathcal{J}_{1}}=\begin{bmatrix}J_{11}&J_{14}\\\
J_{41}&J_{44}\\\ \end{bmatrix}.$

Now

$\displaystyle{\rm tr\ }\overline{\mathcal{J}_{1}}$
$\displaystyle=J_{11}+J_{44}=-(\sigma+\mu)-(t^{\prime}+\mu)<0$
$\displaystyle\det\overline{\mathcal{J}_{1}}$
$\displaystyle=(\sigma+\mu)(t^{\prime}+\mu)-\sigma(1-\lambda)t^{\prime}$
$\displaystyle=\mu(\sigma+\mu+t^{\prime})+\sigma\lambda t^{\prime}>0.$

Thus the eigenvalues of $\overline{\mathcal{J}_{1}}$ have negative real parts.
The remaining three eigenvalues of the Jacobian are the eigenvalues of the
matrix

$\displaystyle\overline{\mathcal{J}_{2}}=\begin{bmatrix}J_{22}&J_{23}&J_{26}\\\
J_{32}&J_{33}&0\\\ 0&J_{63}&J_{66}\end{bmatrix},$

leading to characteristic polynomial

$\displaystyle(\ell-J_{22})(\ell-J_{33})(\ell-J_{66})-J_{23}J_{32}(\ell-
J_{66})-J_{26}J_{32}J_{63}=0,$

which can be rewritten as

$\displaystyle\ell^{3}+\mathcal{A}_{1}\ell^{2}+\mathcal{A}_{2}\ell+\mathcal{A}_{3}=0,$
(14)

where

$\displaystyle\mathcal{A}_{1}$ $\displaystyle=-J_{22}-J_{33}-J_{66},$
$\displaystyle\mathcal{A}_{2}$
$\displaystyle=J_{22}J_{33}-J_{23}J_{32}+J_{66}(J_{22}+J_{33}),$
$\displaystyle\mathcal{A}_{3}$
$\displaystyle=-J_{22}J_{33}J_{66}\left(1-\frac{J_{32}\left(J_{23}J_{66}-J_{63}J_{26}\right)}{J_{22}J_{33}J_{66}}\right).$
(15)

The Routh-Hurwitz criterion says that all roots of the characteristic equation
(14) have negative real parts if and only if $\mathcal{A}_{1}>0$,
$\mathcal{A}_{1}\mathcal{A}_{2}>\mathcal{A}_{3}$ and $\mathcal{A}_{3}>0$ and
if these conditions hold the disease-free steady state is asymptotically
stable. It is clear that the first condition holds but it is not so easy to
see when the second and third conditions hold. It will later be proved
indirectly using the next generation matrix that they hold in this model for
all values of the parameters.

3.1. The next generation matrix

In this section, following the ideas presented in [12], the basic reproduction
ratio for (1)-(6) is derived using the next generation matrix method. We use
the notation of [12]. To apply this method we must choose which of the
unknowns represent groups of infected individuals and which terms in the
equations represent new infections. In fact we choose $E$, $I$ and $C$ to be
the infected variables and the terms which are non-negative and non-linear in
the unknowns to represent new infections. The conditions (A1)-(A5) of [12] are
satisfied. Most of these are rather obvious for this model. The only exception
is (A5) which holds because the quantities corresponding to $J_{23}$ and
$J_{26}$ are zero in the case that new infections are turned off. The matrix
${\mathcal{F}}$ associated with new infections and the matrix ${\mathcal{V}}$
containing the remaining expressions are given by

$\displaystyle{\mathcal{F}}$ $\displaystyle=\begin{bmatrix}0\\\ \beta
SI+\epsilon VI+(\alpha_{1}S+\alpha_{2}V)g(C,\kappa)\\\ 0\\\ 0\\\ 0\\\ 0\\\
\end{bmatrix},$

and

$\displaystyle-{\mathcal{V}}$ $\displaystyle=\begin{bmatrix}0\\\ \xi E+\mu
E\\\ -\xi E+\delta I+dI+\mu I\\\ 0\\\ 0\\\ -\varphi I+\omega C\\\
\end{bmatrix},$

${\mathcal{V}}_{+}$ and ${\\\ mathcal{V}}_{-}$ are the positive and negative
parts of $\mathcal{V}$, respectively. Hence the matrices $F$ and $V$ of [12]
are given by

$\displaystyle F\bigr{\rvert}_{E^{*}}$ $\displaystyle=\begin{bmatrix}0&\beta
S^{*}+\epsilon
V^{*}&\delta_{1n}\left(\alpha_{1}S^{*}+\alpha_{2}V^{*}\right)/\kappa\\\
0&0&0\\\ 0&0&0\\\ \end{bmatrix},\
V\bigr{\rvert}_{E^{*}}=\begin{bmatrix}\xi+\mu&0&0\\\ -\xi&\delta+d+\mu&0\\\
0&-\varphi&\omega\\\ \end{bmatrix},$

The reproduction ratio $R_{0}$ is defined (cf. [6], [12]) to be the spectral
radius of the matrix given by

$\displaystyle
FV^{-1}=\begin{bmatrix}FV^{-1}_{11}&FV^{-1}_{12}&FV^{-1}_{13}\\\ 0&0&0\\\
0&0&0\\\ \end{bmatrix},$

where

$\displaystyle FV^{-1}_{11}$
$\displaystyle=\frac{\xi}{(\xi+\mu)(\mu+\delta+d)}\left[\left(\beta+\frac{\delta_{1n}\alpha_{1}\varphi}{\omega\kappa}\right)S^{*}+\left(\epsilon+\frac{\delta_{1n}\alpha_{2}\varphi}{\omega\kappa}\right)V^{*}\right],$
$\displaystyle FV^{-1}_{12}$
$\displaystyle=\frac{1}{(\mu+\delta+d)}\left[\left(\beta+\frac{\delta_{1n}\alpha_{1}\varphi}{\omega\kappa}\right)S^{*}+\left(\epsilon+\frac{\delta_{1n}\alpha_{2}\varphi}{\omega\kappa}\right)V^{*}\right],$
$\displaystyle FV^{-1}_{13}$
$\displaystyle=\frac{\delta_{1n}}{\omega\kappa}\left(\alpha_{1}S^{*}+\alpha_{2}V^{*}\right)$

The characteristic equation of this matrix is given by

$\displaystyle\det(FV^{-1}-\Sigma I)=0.$

Its roots are the eigenvalues:

$\displaystyle\Sigma_{1}$
$\displaystyle=\frac{\xi}{(\xi+\mu)(\mu+\delta+d)}\left[\left(\beta+\frac{\delta_{1n}\alpha_{1}\varphi}{\omega\kappa}\right)S^{*}+\left(\epsilon+\frac{\delta_{1n}\alpha_{2}\varphi}{\omega\kappa}\right)V^{*}\right],$
(16) $\displaystyle\Sigma_{2}$ $\displaystyle=0,$ (17)
$\displaystyle\Sigma_{3}$ $\displaystyle=0.$ (18)

Thus the basic reproduction ratio, which is associated with the dominant
eigenvalue $\Sigma_{1}$, is

$\displaystyle R_{0}$ $\displaystyle=\frac{\xi
S^{*}}{(\xi+\mu)(\mu+\delta+d)}\left[\left(\beta+\frac{\delta_{1n}\alpha_{1}\varphi}{\omega\kappa}\right)+\left(\epsilon+\frac{\delta_{1n}\alpha_{2}\varphi}{\omega\kappa}\right)\frac{\sigma}{t^{\prime}+\mu}\right],$
(19)
$\displaystyle=\frac{\xi}{\omega(\xi+\mu)(\mu+\delta+d)}\frac{\Lambda(t^{\prime}+\mu)}{\mu(\sigma+t^{\prime}+\mu)+\lambda
t^{\prime}\sigma}\left[\omega\beta+\frac{\delta_{1n}\alpha_{1}\varphi}{\kappa}+\frac{\sigma}{t^{\prime}+\mu}\left(\epsilon\omega+\frac{\delta_{1n}\alpha_{2}\varphi}{\kappa}\right)\right].$

It follows from Theorem 2 of [12] that the disease-free steady state is
asymptotically stable in the case $R_{0}<1$ and unstable if $R_{0}>1$. In fact
looking at the proof reveals that under the assumptions of that theorem the
following stronger statements hold. In the case $R_{0}<1$ all eigenvalues of
the linearization at the disease-free steady state have negative real parts
and in the case $R_{0}>1$ the linearization has an eigenvalue with positive
real part. This gives an indirect proof that the inequalities
$\mathcal{A}_{1}\mathcal{A}_{2}>\mathcal{A}_{3}$ and $\mathcal{A}_{3}>0$ of
the last section hold.

In order to understand the effects of vaccination it is useful to write the
basic reproductive ratio schematically in the form
$R_{0}=A\left(\frac{B\sigma+C}{D\sigma+E}\right)$ where

$\displaystyle A=\frac{\Lambda\xi}{\omega(\xi+\mu)(\mu+\delta+d)},$ (20)
$\displaystyle B=\epsilon\omega+\frac{\delta_{1n}\alpha_{2}\varphi}{\kappa},$
(21) $\displaystyle
C=(\mu+t^{\prime})\left(\omega\beta+\frac{\delta_{1n}\alpha_{1}\varphi}{\kappa}\right),$
(22) $\displaystyle D=\mu+\lambda t^{\prime},$ (23) $\displaystyle
E=\mu(\mu+t^{\prime}).$ (24)

The sign of the derivative of $R_{0}$ with respect to $\sigma$ is equal to
that of $BE-CD$. This last quantity is equal to

$\omega(\mu+t^{\prime})[(\epsilon-\beta)\mu-\lambda
t^{\prime}\beta]+\frac{\delta_{1n}\varphi}{\kappa}(\mu+t^{\prime})[(\alpha_{2}-\alpha_{1})\mu-\lambda
t^{\prime})].$ (25)

Under the assumptions made on the parameters it is negative and so we see that
increasing the vaccination rate decreases $R_{0}$, generalizing a result of
[12].

3.2. Backward bifurcation analysis

The concept of a backward bifurcation is used in the literature on
epidemiological models. It is defined in situations where a definition of the
basic reproduction ratio $R_{0}$ is available. In many models endemic steady
states only exist in the case $R_{0}>1$. We think of the direction of
increasing $R_{0}$ as the forward direction and that is where endemic steady
states occur. There are, however, models where it happens that near $R_{0}=1$
there are endemic steady states with $R_{0}<1$, i.e. in the backward
direction. The steady state bifurcates from the disease-free steady state as
the parameter $R_{0}$ is varied. A commonly occurring case is where there is a
generic transcritical bifurcation for $R_{0}=1$ and this covers both the
forward and backward cases, these being distinguished by the sign of a
parameter $a$. We call this case a generic backward bifurcation. In our model
the qualitative behaviour depends on a parameter which is a natural number
$n$. For $n\neq 2$ we show that any generic transcritical bifurcation must be
a forward bifurcation. Thus a generic backward bifurcation is impossible in
that case. For $n=2$ we show that generic backward bifurcations do occur for
some values of the parameters.

Consider the disease-free steady state $E^{*}=(S^{*},0,0,V^{*},R^{*},0)$ for
the system (1)-(6). We choose $\beta$ as the bifurcation parameter and denote
its value at the bifurcation point where $R_{0}=1$ by $\beta^{*}$. Then

$\displaystyle\beta^{*}=\frac{(\xi+\mu)(\mu+\delta+d)}{\xi
S^{*}}-\frac{\delta_{1n}\alpha_{1}\varphi}{\omega\kappa}-\left(\epsilon+\frac{\delta_{1n}\alpha_{2}\varphi}{\omega\kappa}\right)\frac{\sigma}{t^{\prime}+\mu},$
(26)

where
$S^{*}=\dfrac{\Lambda(t^{\prime}+\mu)}{\mu(\sigma+t^{\prime}+\mu)+\lambda
t^{\prime}\sigma}.$ Note that for fixed values of the other parameters there
is only a choice of $\beta$ for which $R_{0}=1$ if the right hand side of (26)
is positive.

Following the ideas presented by [12], according to center manifold theory, it
is necessary to compute right and left eigenvectors of the Jacobian matrix
evaluated at the disease-free steady state $E^{*}$ and $\beta=\beta^{*}$.
Consider a right eigenvector of the form
$\underline{w}=(w_{1},w_{2},w_{3},w_{4},w_{5},w_{6})^{T}$. Thus the system
leads to

$\displaystyle-(\sigma+\mu)w_{1}-\beta^{*}S^{*}w_{3}+(1-\lambda)t^{\prime}w_{4}-\frac{\alpha_{1}\delta_{1n}}{\kappa}S^{*}w_{6}$
$\displaystyle=0,$ (27) $\displaystyle-(\xi+\mu)w_{2}+(\beta^{*}S^{*}+\epsilon
V^{*})w_{3}+\frac{\delta_{1n}}{\kappa}(\alpha_{1}S^{*}+\alpha_{2}V^{*})w_{6}$
$\displaystyle=0,$ (28) $\displaystyle\xi w_{2}-(\delta+d+\mu)w_{3}$
$\displaystyle=0,$ (29) $\displaystyle\sigma w_{1}-\epsilon
V^{*}w_{3}-(t^{\prime}+\mu)w_{4}-\frac{\delta_{1n}\alpha_{2}}{\kappa}V^{*}w_{6}$
$\displaystyle=0,$ (30) $\displaystyle\delta\omega_{3}+\lambda
t^{\prime}w_{4}-\mu w_{5}$ $\displaystyle=0,$ (31) $\displaystyle\varphi
w_{3}-\omega w_{6}$ $\displaystyle=0.$ (32)

Using Eq. (29) and (32), we obtain

$\displaystyle w_{3}=\dfrac{\xi w_{2}}{\delta+d+\mu}\quad\textrm{and}\quad
w_{6}=\dfrac{\varphi\xi}{\omega(\delta+d+\mu)}w_{2}.$ (33)

Using these in Eq. (28), we find

$\displaystyle\left[\beta^{*}-\frac{(\delta+d+\mu)(\xi+\mu)}{\xi
S^{*}}+\frac{\delta_{1n}\alpha_{1}\varphi}{\omega\kappa}+\left(\epsilon+\frac{\delta_{1n}\alpha_{2}\varphi}{\omega\kappa}\right)\frac{\sigma}{t^{\prime}+\mu}\right]w_{2}=0,$
(34)

Thus we see that (26) is a necessary condition for there to be a vector in the
kernel with $w_{2}\neq 0$. Note that if $w_{2}>0$ then $w_{3}$ and $w_{6}$ are
positive, as they must be as a consequence of the general theory. Besides,
using Eqs. (27) and (30), one obtains

$\displaystyle-(\sigma+\mu)w_{1}+(1-\lambda)t^{\prime}w_{4}$
$\displaystyle=\left(\beta^{*}+\frac{\delta_{1n}\alpha_{1}\varphi}{\omega\kappa}\right)\frac{\xi
S^{*}}{(\delta+d+\mu)}w_{2},$ (35) $\displaystyle\sigma
w_{1}-(t^{\prime}+\mu)w_{4}$
$\displaystyle=\left(\epsilon+\frac{\delta_{1n}\alpha_{2}\varphi}{\omega\kappa}\right)\frac{\xi
V^{*}}{(\delta+d+\mu)}w_{2},$ (36)

respectively. That leads to

$\displaystyle w_{1}=-\dfrac{\dfrac{\xi
w_{2}}{\delta+d+\mu}\left[\left(\beta^{*}+\dfrac{\delta_{1n}\alpha_{1}\varphi}{\omega\kappa}\right)(t^{\prime}+\mu)S^{*}+\left(\epsilon+\dfrac{\delta_{1n}\alpha_{2}\varphi}{\omega\kappa}\right)(1-\lambda)t^{\prime}V^{*}\right]}{\mu(\sigma+t^{\prime}+\mu)+\sigma\lambda
t^{\prime}}<0,$ (37)

for $\lambda\leq 1$ and

$\displaystyle w_{4}=-\dfrac{\dfrac{\xi
w_{2}}{\delta+d+\mu}\left[\left(\beta^{*}+\dfrac{\delta_{1n}\alpha_{1}\varphi}{\omega\kappa}\right)\sigma
S^{*}+\left(\epsilon+\dfrac{\delta_{1n}\alpha_{2}\varphi}{\omega\kappa}\right)(\sigma+\mu)V^{*}\right]}{\mu(\sigma+t^{\prime}+\mu)+\sigma\lambda
t^{\prime}}<0.$ (38)

Furthermore using Eq. (31):

$\displaystyle w_{5}=\dfrac{\xi w_{2}}{\delta+d+\mu}\left[\delta-\lambda
t^{\prime}\dfrac{\left(\beta^{*}+\dfrac{\delta_{1n}\alpha_{1}\varphi}{\omega\kappa}\right)\sigma
S^{*}+\left(\epsilon+\dfrac{\delta_{1n}\alpha_{2}\varphi}{\omega\kappa}\right)(\sigma+\mu)V^{*}}{\mu(\sigma+t^{\prime}+\mu)+\sigma\lambda
t^{\prime}}\right].$ (39)

In a similar manner, a left eigenvector can be written in the form
$\underline{v}=(v_{1},v_{2},v_{3},v_{4},v_{5},v_{6})$ for which

$\displaystyle v_{5}$ $\displaystyle=0,$ (40)
$\displaystyle-(\sigma+\mu)v_{1}+\sigma v_{4}$ $\displaystyle=0,$ (41)
$\displaystyle-(\xi+\mu)v_{2}+\xi v_{3}$ $\displaystyle=0,$ (42)
$\displaystyle-\beta^{*}S^{*}v_{1}+(\beta^{*}S^{*}+\epsilon
V^{*})v_{2}-(\delta+d+\mu)v_{3}-\epsilon V^{*}v_{4}+\varphi v_{6}$
$\displaystyle=0,$ (43)
$\displaystyle(1-\lambda)t^{\prime}v_{1}-(t^{\prime}+\mu)v_{4}$
$\displaystyle=0,$ (44)
$\displaystyle-\frac{\delta_{1n}\alpha_{1}}{\kappa}S^{*}v_{1}+\frac{1}{\kappa}(\delta_{1n}\alpha_{1}S^{*}+\delta_{1n}\alpha_{2}V^{*})v_{2}-\frac{\delta_{1n}\alpha_{2}}{\kappa}V^{*}v_{4}-\omega
v_{6}$ $\displaystyle=0.$ (45)

Using (41) and (44), we find

$\displaystyle\left(\sigma\lambda
t^{\prime}+\mu(\sigma+\mu+t^{\prime})\right)v_{4}=0,$ (46)

leading to $v_{4}=0$ and thus $v_{1}=0$. In addition, from Eqs. (42) and (45),
we obtain

$\displaystyle
v_{3}=\left(1+\frac{\mu}{\xi}\right)v_{2},\quad\textrm{and}\quad
v_{6}=\frac{\delta_{1n}}{\kappa\omega}\left(\alpha_{1}S^{*}+\alpha_{2}V^{*}\right)v_{2}.$
(47)

Thus the left eigenvector becomes

$\displaystyle\underline{v}=\left(0,v_{2},\left(1+\frac{\mu}{\xi}\right)v_{2},0,0,\frac{\delta_{1n}}{\kappa\omega}(\alpha_{1}S^{*}+\alpha_{2}V^{*})v_{2}\right),$
(48)

and we find

$\displaystyle\underline{w}\cdot\underline{v}=\left(1+\frac{\xi+\mu}{\delta+d+\mu}+\frac{\delta_{1n}\varphi\xi}{\kappa\omega^{2}(\delta+d+\mu)}\right)w_{2}v_{2}>0.$
(49)

Let $a$ be the bifurcation coefficient introduced in [12]. Considering the
model (1)-(6) in the form $\dot{x_{i}}=f_{i}(x_{i}),\;\;i=\\{1,2,3,4,5,6\\}$,
it is given by

$\displaystyle
a=\sum\limits_{k,i,j=1}^{6}v_{k}w_{i}w_{j}\frac{\partial^{2}f_{k}}{\partial
x_{i}\partial x_{j}}(0,0).$ (50)

Using Eqs. (48) and (50)

$\displaystyle a=$ $\displaystyle
v_{2}\sum\limits_{i,j=1}^{6}w_{i}w_{j}\frac{\partial^{2}f_{2}}{\partial
x_{i}\partial
x_{j}}+\left(1+\frac{\mu}{\xi}\right)v_{2}\sum\limits_{i,j=1}^{6}w_{i}w_{j}\frac{\partial^{2}f_{3}}{\partial
x_{i}\partial x_{j}}$
$\displaystyle+\frac{\delta_{1n}(\alpha_{1}x_{1}^{*}+\alpha_{2}x_{4}^{*})}{\kappa\omega}v_{2}\sum\limits_{i,j=1}^{6}w_{i}w_{j}\frac{\partial^{2}f_{6}}{\partial
x_{i}\partial x_{j}},$

where

$\displaystyle f_{2}$ $\displaystyle=\beta^{*}x_{1}x_{3}+\epsilon
x_{3}x_{4}-\xi x_{2}-\mu
x_{2}+\alpha_{1}x_{1}g(x_{6},\kappa)+\alpha_{2}x_{4}g(x_{6},\kappa),$
$\displaystyle f_{3}$ $\displaystyle=\xi x_{2}-(\delta+d+\mu)x_{3},$
$\displaystyle f_{6}$ $\displaystyle=\varphi x_{3}-\omega x_{6}.$

Note that $\frac{\partial g}{\partial C}(0,\kappa)=\kappa^{-1}$ for $n=1$ and
zero otherwise. $\frac{\partial^{2}g}{\partial C^{2}}(0,\kappa)$ is equal to
$-2\kappa^{-2}$ for $n=1$, $2\kappa^{-1}$ for $n=2$ and zero otherwise. Since
second derivatives of $f_{3}$ and $f_{6}$ with respect to
$x_{i},\;\;i=\\{1,2,3,4,5,6\\}$ are always zero;

$\displaystyle a=2v_{2}\left[w_{1}w_{3}\beta^{*}+w_{3}w_{4}\epsilon\right.$
$\displaystyle\left.+(\alpha_{1}w_{1}+\alpha_{2}w_{4})w_{6}\frac{\partial
g}{\partial
C}(0,\kappa)+(\alpha_{1}x_{1}^{*}+\alpha_{2}x_{4}^{*})w_{6}^{2}\frac{\partial^{2}g}{\partial
C^{2}}(0,\kappa)\right].$ (51)

We want to determine the sign of $a$ and since $v_{2}>0$ this is the same as
that of the expression in square brackets. The first two summands are
negative. Now $(w_{1}+w_{4})w_{6}<0$ and $w_{6}^{2}>0$. Hence in the case
$n=1$ we see that all summands are negative and hence $a<0$. It follows that
in that case the conditions for a backward bifurcation given in [12] cannot be
satisfied. The same conclusion is obtained in the case $n\geq 3$. In the
exceptional case $n=2$ the last summand is positive and so we investigate
further whether a backward bifurcation can take place in that case. In the
notation of [12] we have

$\displaystyle
b=\sum\limits_{k,i=1}^{6}w_{k}w_{i}\frac{\partial^{2}f_{k}}{\partial
x_{i}\partial\beta}(0,0).$ (52)

Since only $f_{1}$ and $f_{2}$ involve the parameter $\beta$, derivatives with
respect to $x_{1}$ and $x_{3}$ are non-zero. Here

$\displaystyle b=w_{3}x_{1}^{*}(w_{2}-w_{1})>0.$ (53)

It follows from [12] that there is a backward bifurcation precisely when the
right hand side of (26) is positive and

$2(\alpha_{1}S^{*}+\alpha_{2}V^{*})w_{6}^{2}\kappa^{-1}>-w_{3}(w_{1}\beta^{*}+w_{4}\epsilon).$
(54)

Note that in the case $n>1$ the expression for $\beta^{*}$ simplifies to

$\displaystyle\beta^{*}=\frac{(\xi+\mu)(\mu+\delta+d)\left[\mu(\sigma+t^{\prime}+\mu)+\lambda
t^{\prime}\sigma\right]-\Lambda\xi\epsilon\sigma}{\Lambda\xi(t^{\prime}+\mu)}.$
(55)

Do there exist values of the parameters for which these conditions are
satisfied? To investigate this we substitute the expressions for $w_{6}$,
$w_{3}$, $w_{1}$ and $w_{4}$ into (54). The result is

$\displaystyle
2(\alpha_{1}S^{*}+\alpha_{2}V^{*})\dfrac{\varphi^{2}}{\omega^{2}\kappa}$
$\displaystyle>\frac{\Lambda[(\beta^{*})^{2}(t^{\prime}+\mu)^{2}+\beta^{*}\epsilon((2-\lambda)t^{\prime}+\mu)\sigma+\epsilon^{2}t^{\prime}\sigma(\sigma+\mu)]}{[\mu(\sigma+t^{\prime}+\mu)+\sigma\lambda
t^{\prime}]^{2}}.$ (56)

Using the expressions for $S^{*}$ and $V^{*}$ this can be simplified to

$\displaystyle
2(\alpha_{1}(t^{\prime}+\mu)+\alpha_{2}\sigma)\dfrac{\varphi^{2}}{\omega^{2}\kappa}$
$\displaystyle>\frac{[(\beta^{*})^{2}(t^{\prime}+\mu)^{2}+\beta^{*}\epsilon((2-\lambda)t^{\prime}+\mu)\sigma+\epsilon^{2}t^{\prime}\sigma(\sigma+\mu)]}{\mu(\sigma+t^{\prime}+\mu)+\sigma\lambda
t^{\prime}}.$ (57)

Thus it is clear that if $\alpha_{1}$ or $\alpha_{2}$ is made large enough
while the other parameters are kept fixed then there is a backward
bifurcation. It is important to note that the parameters given for the case
$n=2$ in Table 3 satisfy the conditions for a backward bifurcation in (57).
These results with now be summed up.

Theorem 3.1 If $n=2$ and the parameters in the system (1)-(6) satisfy the
inequality 57 with the quantity $\beta^{*}$ defined by 55 being positive then
the parameter $a$ of [12] is positive and a generic backward bifurcation
occurs. There exist parameters for which these conditions are satisfied. If
$n\neq 2$ the condition $a>0$ is never satisfied.

The centre manifold at the bifurcation point is one-dimensional. Since
$v_{3}>0$ we can use $I$ as a parameter on the centre manifold. The
restriction of the dynamical system to the centre manifold is of the form
$\dot{I}=f(I,\beta)$, where we have suppressed the dependence on the
parameters other than $\beta$ in the notation. With this notation the sign of
$a$ is equal to that of $\frac{\partial^{2}f}{\partial I^{2}}$ while that of
$b$ is equal to that of $\frac{\partial^{2}f}{\partial I\partial\beta}$. We
have $f(0,\beta)=0$ for all $\beta$. Moreover there is a curve $c(\beta)$ of
steady states with $c(\beta^{*})=0$. The sign of $c^{\prime}(\beta^{*})$ is
equal to that of $a$.

4\. Endemic steady states

In this section we consider endemic steady states, i.e. those where all the
unknowns are positive. It follows from (4) that at a steady state
$V=\frac{\sigma}{\epsilon I+(t^{\prime}+\mu)+\alpha_{2}g}S$. Substituting this
into (1) gives

$\displaystyle S=\frac{\Lambda}{\beta
I+\sigma+\alpha_{1}g+\mu}\left[1-\frac{(1-\lambda)t^{\prime}\sigma}{(\beta
I+\sigma+\alpha_{1}g+\mu)(\epsilon
I+(t^{\prime}+\mu)+\alpha_{2}g)}\right]^{-1}$
$\displaystyle=\frac{\Lambda(\epsilon
I+(t^{\prime}+\mu)+\alpha_{2}g(C))}{(\beta
I+\sigma+\alpha_{1}g(C)+\mu)(\epsilon
I+(t^{\prime}+\mu)+\alpha_{2}g(C))-(1-\lambda)t^{\prime}\sigma}$ (58)

Note also that due to (6) we have $C=\frac{\varphi}{\omega}I$ and that due to
(3) we have $E=\frac{\xi}{\delta+d+\mu}I$. These relations allow $S$, $E$, $V$
and $C$ to be expressed in terms of $I$. Thus substituting them into (2) gives
an equation for $I$ alone. Each summand contains a factor $I$ and for an
endemic steady state this can be cancelled. There remains

$0=\beta S+\epsilon
V-\frac{(\xi+\mu)(\delta+d+\mu)}{\xi}+(\alpha_{1}S+\alpha_{2}V)\frac{(\varphi/\omega)^{n}I^{n-1}}{(\varphi/\omega)^{n}I^{n}+\kappa}$
(59)

Denote the expression in the denominator of (Dynamics of a mathematical model
of virus spreading incorporating the effect of a vaccine) by $Z$. Multiplying
(59) by $Z(\epsilon
I+(t^{\prime}+\mu)+\alpha_{2}g)[(\varphi/\omega)^{n}I^{n}+\kappa]$ gives a
polynomial equation $p(I)=0$ for $I$. Endemic steady states are in one to one
correspondence with positive roots of $p$.

Since we are most interested in the case with backward bifurcations we now
restrict to the case $n=2$, where there are considerable simplifications in
these formulae. It is clear that $p(I)\to-\infty$ as $I\to\infty$. Moreover
$p(0)$ is equal to a positive factor times

$\beta S^{*}+\epsilon V^{*}-\frac{(\xi+\mu)(\delta+d+\mu)}{\xi}.$ (60)

We see that the sign of $p(0)$ is the same as that of $R_{0}-1$. As $\beta$
increases through $\beta^{*}$ the sign of $p(0)$ changes from negative to
positive. When $p^{\prime}(0)\neq 0$ we have a generic transcritical
bifurcation (cf. [28], Section 3.1). The sign of $p^{\prime}(0)$ is the same
as that of $c^{\prime}(\beta^{*})$ in the discussion of the centre manifold in
the previous section. If $p^{\prime}(0)<0$ then for $\beta$ slightly greater
than $\beta^{*}$ there exists a positive root of $p$ close to zero and hence a
positive steady state close to the disease free steady state. This corresponds
to a forward bifurcation. If, on the other hand, $p^{\prime}(0)>0$ then for
$\beta$ slightly less than $\beta^{*}$ there exists a positive root of $p$
close to zero and hence a positive steady state close to the disease free
steady state. This corresponds to a backward bifurcation. In this case $p$ is
positive for $I$ slightly larger than its value $I_{1}$ at that steady state.
By the intermediate value theorem there must exist some $I_{2}>I_{1}$ with
$p(I_{2})=0$ and hence a second positive steady state. The direction of the
flow on the centre manifold shows that the steady state with $I=I_{1}$ is
unstable. The stability of the steady state with $I=I_{2}$ cannot be
determined by the arguments we have presented.

5\. Simulations

In a broad context, the process of mathematical modelling and data fitting
revolves around formulating mathematical models that describe real-world
phenomena and then adjusting the parameters of these models to best match
observed data. Therefore, the aim is to capture the underlying relationships
and behavior of the system being studied and use the available data to
validate the model. In our study, we ensure the validation of the mathematical
model of the COVID-19 outbreak by using the data fitting of the model
regarding the observed data. However, the available data is scarce (only the
actual data for infectious and vaccinated people are almost certain).
Vaccination is globally considered to be the most effective solution for
infectious diseases such as the recent COVID-19 outbreak. In this section, as
an example, model results and observed data for the vaccinated class are
compared for COVID-19 scenarios in Turkey. Some of the realistic parameters
are taken from the literature, see the related references in Table 3 for a
detailed discussion of parameter choices. The remaining seven parameters are
estimated, by fitting the vaccinated compartment generated from the system
(1)-(6) to the observed number of COVID-19 vaccinated individuals using
standard model-fitting procedures.

The least squares method is the process of finding the best-fitting curve or
line of best fit for a set of data points by reducing the sum of the squares
of the offsets (residual part) of the points from the curve. The vector
consisting of seven parameters
$p=(\beta,\epsilon,\alpha_{1},\alpha_{2},\xi,\delta,\sigma)$ can be estimated
via parameter estimation. In this context, the model given by (1)-(6) is
evaluated by considering a non-linear least squares problem with positive
constraints, where the best fitting curve can be found for a small data set of
vaccinated class by minimising the sum of squares of the deviations of data
points from the plotted graph [6, 8]. This may be described as

$S(p)=\sum\left(V_{i}-\mathcal{F}(x_{i},p)\right)^{2},$ (61)

where $V_{i}$ represents the data set and $\mathcal{F}(x_{i},p)$ denotes the
model result with for a vector of unknowns $p$. To minimise the function
$S(p)$, the non-linear least square minimization routine lsqcurvefit of MATLAB
is used [17]. Parameters obtained from this approach are given in Table 3.
Besides MATLAB’s standard ode45 solver [17] is applied for numerical
integration of the system (1)-(6) with suitable initial conditions provided in
the text.

Numerical simulations of the model (1)-(6) can be shown with the parameters
given in Table 3 for $n=1$. An example data set of vaccinated people during
the COVID$-19$ outbreak in Turkey is taken from the World Health Organisation.
In Fig. 2, the results are shown for model (1)-(6) with $n=1$ fitted to the
data of individuals vaccinated between June $10$, $2021$ and August $8$,
$2021$. The total population is taken $83$ million [1].

h

Figure 2. Vaccination component for the model for the case with $n=1$ is
compared with real data for the period of June $10$, $2021$ and August $8$,
$2021$. Parameters are given in Table 3. Initial conditions are chosen as
$S_{0}=61098000,V_{0}=18500000,E_{0}=2200000,R_{0}=2000,C_{0}=20000$. The
black line denotes the model result and the red stars represent daily
vaccinated cases.

In Fig. 3, success of parameter estimation is again demonstrated. Here
numerical simulations of the model (1)-(6) with $n=2$ is shown with the
parameters given in Table 3 and the resulting outcome is compared with the
data set of vaccinated people between June $10$, $2021$ and August $8$,
$2021$. As seen, in Figs. 2 and 3, the black curve corresponding to the model
becomes flattened after around day $40$ and this agrees with the real data,
where the real number of vaccinated people rapidly increases from $18967237$
to $41726338$ and it is smoothed roughly about July $20$.

Figure 3. Vaccination component for the model for the case with $n=2$ is compared with real data for the period of June $10$, $2021$ and August $8$, $2021$. Parameters are given in Table 3. Initial conditions are chosen as $S_{0}=61098000,V_{0}=18500000,E_{0}=2200000,R_{0}=2000,C_{0}=20000$. The black line denotes the model result and the red stars represent daily vaccinated cases.  Table 3. Estimated parameters of the model. Parameters | Value ($n=1$) | Value ($n=2$) | Unit | Source
---|---|---|---|---
$\Lambda$ | $3032$ | $3032$ | day-1 | assumed based on [5, 14, 15]
$\beta$ | $1.0257\times 10^{-8}$ | $1.004\times 10^{-8}$ | day-1 | estimated
$\mu$ | $3.653\times 10^{-5}$ | $3.653\times 10^{-5}$ | day-1 | assumed based on [1, 6, 15]
$\epsilon$ | $1\times 10^{-8}$ | $1\times 10^{-8}$ | day-1 | estimated
$t^{\prime}$ | $0.0055$ | $0.0055$ | day-1 | assumed based on [3]
$\lambda$ | $0.8$ | $0.8$ | day-1 | assumed based on [3, 4, 21]
$d$ | $0.1$ | $0.1$ | day-1 | assumed based on [6, 23]
$\alpha_{1}$ | $0.00041$ | $0.0001$ | day-1 | estimated
$\alpha_{2}$ | $0.00031$ | $0.0001$ | day-1 | estimated
$\omega$ | $5$ | $5$ | day-1 | [6]
$\varphi$ | $2$ | $2$ | day-1 | [6]
$\kappa$ | $20000$ | $20000$ | copies $/$day | assumed based on [13, 6]
$\xi$ | $0.01004$ | $0.01133$ | day-1 | estimated
$\delta$ | $0.19999$ | $0.2$ | day-1 | estimated
$\sigma$ | $0.02136$ | $0.02126$ | day-1 | estimated

6\. Sensitivity analysis

Since varying the parameters may have a significant impact on the model
output, one can perform a sensitivity analysis of the dynamical model to
determine which model parameters are more influential for the dynamics. The
parameters associated with the basic reproduction ratio $R_{0}$ have
particular importance for the robustness of the model. In this context, the
aim of the sensitivity analysis is to identify the most substantial parameter
in the model for disease transmission.

Following the ideas presented in [7, 18], the sensitivity analysis can be
performed based on the basic reproduction ratio as

$\mathcal{S}_{i}^{\mathcal{P}}=\frac{\mathcal{P}}{R_{0}}\frac{\partial
R_{0}}{\partial\mathcal{P}},$ (62)

where $\mathcal{P}$ represents a generic parameter in the model (1)-(6).

The sensitivity indices of the system parameters in Table 3 are demonstrated
in 4, and also in Figs. 4(a) for $n=1$ and 4(b) for $n=2$.

As seen from Table 4 and Fig. 4, the model (1)-(6) is highly sensitive to the
parameters $\lambda$, $t^{\prime}$ and $\delta$. Thus it can be concluded that
an increase in these parameters diminishes the basic reproduction ratio
$R_{0}$ for both $n=1$ and $n=2$. The significance of some parameters may be
different between Case $1$ and Case $2$. For example, although the increase in
the parameter $\epsilon$, i.e. the rate at which a vaccinated individual
becomes exposed after being in contact with an infected individual, leads to
an essential stimulus for the basic reproduction ratio for Case $2$, yet it
has a much smaller impact on the $R_{0}$ for the case $n=1$.

Table 4. Sensitivity indices of basic reproduction ratio of model (1)-(6), considering at the baseline parameters provided in Table 3. number for $n=1$ and $n=2$ Parameters | Index ($n=1$) | Index ($n=2$)
---|---|---
$\beta$ | $0.1276764$ | $0.2089141$
$\epsilon$ | $0.4754643$ | $0.7910858$
$t\prime$ | $-0.7636042$ | $-0.7843018$
$\lambda$ | $-0.9897355$ | $-0.9897256$
$d$ | $-0.3333038$ | $-0.3332927$
$\alpha_{1}$ | $0.1020714$ | $0$
$\alpha_{2}$ | $0.2947878$ | $0$
$\xi$ | $0.0036252$ | $0.0032138$
$\delta$ | $-0.6665743$ | $-0.6665855$
$\sigma$ | $-0.2276182$ | $-0.2067744$

---
(a)

---
(b)

Figure 4. Sensitivity indices based on basic reproduction ratio $R_{0}$ with
respect to various parameters given for Cases $n=1$ (a) and $n=2$ (b) in the
model (1)-(6) for Turkey. Parameter values given in Table 3 are considered for
both cases.

Since the aim of our work is to further broaden the current knowledge of the
modelling of recent COVID-19 outbreak with vaccination, we focus on the role
of two parameters regarding vaccination in the model. In Fig. 5, time
simulations of the Susceptible compartment over a period of $180$ days are
presented for various rates of $\lambda=0.2,0.8,1,4$ (a) and
$\sigma=0.01,0.02136,0.09$ (b), respectively denotes the efficiency of the
vaccine and the vaccination rate of susceptible (after the first shot). As is
seen from Fig. 5 (a,b), with the increase in $\lambda$ and $\sigma$,
susceptible individuals diminish at an earlier time and enter the vaccinated
class.

---
(a)

---
(b)

Figure 5. Plots for Susceptible population compartment for different values of
vaccine efficiency $\lambda=0.2,0.8,1.4$ and for different values of
vaccination rate of susceptible individuals after first shot,
$\sigma=0.01,0.02136,0.09$ with $n=1$.

Figure 6 demonstrates time simulations of the recovered class over a period of
$360$ days for various rates of $\lambda=0.2,0.8,1,4$ and
$\sigma=0.01,0.02136,0.09$.

---
(a)

---
(b)

Figure 6. Plots for Recovery compartment for different values of vaccine
efficiency $\lambda=0.2,0.8,1.4$; and for different values of vaccination rate
of susceptible individuals after first shot, $\sigma=0.01,0.02136,0.09$ with
$n=1$.

Conclusion and outlook

In this paper a model for an epidemic with a partially effective vaccination
and infection by virus in the environment is studied. Two different
implementations of the idea of a partially effective vaccination are included
and which of these is chosen does not seem to have an essential influence on
the qualitative behaviour of the solutions. In modelling infections coming
from the environment we used a phenomenological model for the force of
infection containing an integer $n\geq 1$ as a parameter. We discovered that
the choice $n=2$ leads to the occurrence of backward bifurcations while other
choices of $n$ do not. Thus here there is a major qualitative difference. On
the other hand, fitting to real data for COVID-19 in Turkey showed that both
the cases $n=1$ and $n=2$ worked and there was no clear indication that one of
these choices was better than the other according to that criterion.

In the case $n=2$ it was shown that for certain values of the parameters there
exists an endemic steady state although $R_{0}<1$. It was also shown that in
this situation there exist more than one endemic steady state. At least one of
the positive steady states is unstable. This confirms rigorously that certain
aspects of the usual picture of a backward bifurcation are present in this
model. An aspect of this picture which was not reproduced here is that the
other positive steady state should be stable. It would be desirable to prove a
stability statement of this kind analytically.

An outstanding challenge is to provide a mechanistic derivation of the
response function for infections coming from the environment. If this could be
done then it would help to decide which value of $n$ in the Ansatz we used is
more appropriate for modelling a given disease or whether, indeed, a different
type of function would give better results.

††footnotetext: Funding: Aytül Gökçe was partially supported by a grant from
the Niels Hendrik Abel Board.

## References

  * [1] Türkiye İstatistik Kurumu (TUİK), https://www.tuik.gov.tr/
  * [2] Anderson, G., Reiter, R.J.: Melatonin: Roles in influenza, COVID-19, and other viral infections. Reviews in Medical Virology 30, e2019 (2020)
  * [3] Angeli, M., Neofotistos, G., Mattheakis, M., Kaxiras, E.: Modeling the effect of the vaccination campaign on the COVID-19 pandemic. Chaos, Solitons & Fractals 154, 111621 (2022)
  * [4] Baden, L.R., El Sahly, H.M., Essink, B., Kotloff, K., Frey, S., Novak, R., Diemert, D., Spector, S.A., Rouphael, N., Creech, C.B., et al.: Efficacy and safety of the mRNA-1273 SARS-CoV-2 vaccine. New England Journal of Medicine 384(5), 403–416 (2021)
  * [5] Bajiya, V.P., Bugalia, S., Tripathi, J.P.: Mathematical modeling of COVID-19: Impact of non-pharmaceutical interventions in India. Chaos: An Interdisciplinary Journal of Nonlinear Science 30(11), 113143 (2020)
  * [6] Bulut, H., Gölgeli, M., Atay, F.M.: Modelling personal cautiousness during the COVID-19 pandemic: a case study for Turkey and Italy. Nonlinear Dynamics 105(1), 957–969 (2021)
  * [7] Chitnis, N., Hyman, J.M., Cushing, J.M.: Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model. Bulletin of Mathematical Biology 70, 1272–1296 (2008)
  * [8] Coleman, T.F., Li, Y.: An interior trust region approach for nonlinear minimization subject to bounds. SIAM Journal on Optimization 6(2), 418–445 (1996)
  * [9] Dawes, J.H.P., Souza, M.: A derivation of Holling’s type I, II and III functional responses in predator–prey systems. Journal of Theoretical Biology 327, 11–22 (2013)
  * [10] Diekmann, O., Heesterbeek, J.A.P.: Mathematical Epidemiology of Infectious Diseases. Wiley (2000)
  * [11] Dietz, K.: Overall population patterns in the transmission cycle of infectious disease cycles, in Population Biology of Infectious Diseases (Eds. R. M. Anderson and R. M. May). Springer (1982)
  * [12] Van den Driessche, P.,Watmough, J.: Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical Biosciences 180(1-2), 29–48 (2002)
  * [13] Drosten, C., Chiu, L.L., Panning, M., Leong, H.N., Preiser, W., Tam, J.S., G¨unther, S., Kramme, S., Emmerich, P., Ng, W.L., et al.: Evaluation of advanced reverse transcription-PCR assays and an alternative PCR target region for detection of severe acute respiratory syndrome-associated coronavirus. Journal of Clinical Microbiology 42(5), 2043–2047 (2004)
  * [14] Garba, S.M., Safi, M.A., Gumel, A.B.: Cross-immunity-induced backward bifurcation for a model of transmission dynamics of two strains of influenza. Nonlinear Analysis: Real World Applications 14(3), 1384–1403 (2013)
  * [15] Gumel, A.B., McCluskey, C.C., Watmough, J.: (2006). An SVEIR model for assessing potential impact of an imperfect anti-SARS vaccine. Mathematical Biosciences and Engineering 105(1), 957–969 (2021)
  * [16] Holling, C.S.: Some characteristics of simple types of predation and parasitism. The Canadian Entomologist 91(7), 385–398 (1959)
  * [17] The MathWorks Inc.: MATLAB version: 9.13.0 (R2020b) (2020). https://www.mathworks.com
  * [18] Martcheva, M.: An Introduction to Mathematical Epidemiology. Springer (2015)
  * [19] Martcheva, M.: Methods for deriving necessary and sufficient conditions for backward bifurcation. Journal of Biological Dynamics 13, 538–566 (2019)
  * [20] Park, S.W., Cornforth, D.M., Dushoff, J., Weitz, J.S.: The time scale of asymptomatic transmission affects estimates of epidemic potential in the COVID-19 outbreak. Epidemics 31, 100392 (2020)
  * [21] Polack, F.P., Thomas, S.J., Kitchin, N., Absalon, J., Gurtman, A., Lockhart, S., Perez, J.L., P´erez Marc, G., Moreira, E.D., Zerbini, C., et al.: Safety and efficacy of the BNT162b2 mRNA COVID-19 vaccine. New England Journal of Medicine 383(27), 2603–2615 (2020)
  * [22] Rendall, A.D.: Mathematics of the NFAT signalling pathway. SIAM Journal on Applied Dynamical Systems 11, 988–1006 (2012)
  * [23] Ritchie, H., Ortiz-Ospina, E., Beltekian, D., Mathieu, E., Hasell, J., Macdonald, B., Giattino, C., Roser, M., Breck Yunits, A., Gavrilov, D., et al.: Mortality risk of COVID-19. Our World in Data [Internet].[cited 5 May 2020]. Available: https://ourworldindata.org/mortality-risk-covid#the-case-fatality-rate (2020)
  * [24] Ross, R.: An application of the theory of probability to the theory of a priori pathometry. Part I. Proceedings of the Royal Society of London. Series A 92, 204–230 (1916)
  * [25] Ross, R., Hudson, H.P.: An application of the theory of probability to the theory of a priori pathometry. Part II. Proceedings of the Royal Society of London. Series A 93, 212–225 (1917)
  * [26] Ross, R., Hudson, H.P.: An application of the theory of probability to the theory of a priori pathometry. Part III. Proceedings of the Royal Society of London. Series A 93, 225–240 (1917)
  * [27] Sunyok, C.J., Fox, L., Ritchie, H., Lanzas, C., Lenhart, S.: Mathematically modelling the effect of touch frequency on the environmental transmission of Clostridioides difficile in healthcare settings. Mathematical Biosciences 340, 108666 (2021)
  * [28] Wiggins, S.: Introduction to Applied Nonlinear Dynamical Systems and Chaos. Springer (1990)
  * [29] Zhao, J., Eisenberg, J.E., Spicknall, I.H., Li, S., Koopman, J.S.: Model analysis of fomite mediated influenza transmission. PloS ONE 7, e51984 (2012)

††footnotetext: ∗(Corresponding author)