1cm1cm1cm1cm

Mutations make pandemics worse or better: modeling SARS-CoV-2 variants and
imperfect vaccination

Sarita Bugalia1, Jai Prakash Tripathi1, Hao Wang2

1 Department of Mathematics, Central University of Rajasthan,
Bandar Sindri, Kishangarh-305817, Ajmer, Rajasthan, India
2 Department of Mathematical and Statistical Sciences,
University of Alberta, Edmonton AB T6G 2G1, Canada

###### Abstract

COVID-19 is a respiratory disease triggered by an RNA virus inclined to
mutations. Since December 2020, variants of COVID-19 (especially Delta and
Omicron) continuously appeared with different characteristics that influenced
death and transmissibility emerged around the world. To address the novel
dynamics of the disease, we propose and analyze a dynamical model of two
strains, namely native and mutant, transmission dynamics with mutation and
imperfect vaccination. It is also assumed that the recuperated individuals
from the native strain can be infected with mutant strain through the direct
contact with individual or contaminated surfaces or aerosols. We compute the
basic reproduction number for each strain independently and take the maximum
for $R_{0}$. We prove the nonexistence of backward bifurcation using the
center manifold theory, and global stability of disease-free equilibrium when
the basic reproduction number $R_{0}<1$, that is, vaccine is effective to
eliminate the native and mutant strains even if it cannot provide full
protection. Hopf bifurcation appears when the endemic equilibrium loses its
stability. An intermediate mutation rate $\nu_{1}$ leads to oscillations. When
$\nu_{1}$ increases over a threshold, the system regains its stability and
exhibits an interesting dynamics called endemic bubble. An analytical
expression for vaccine-induced herd immunity is derived. The epidemiological
implication of the herd immunity threshold is that the disease can be
effectively eradicated if the minimum herd immunity threshold is attained in
the community. Furthermore, the model is parameterized using the Indian data
of the cumulative number of confirmed cases and deaths of COVID-19 from March
1 to September 27 in 2021, using MCMC method. The cumulative cases and deaths
can be reduced by increasing the vaccine efficacies to both native and mutant
strains. We observe that by considering the vaccine efficacy to native strain
as 90%, the cumulative cases and deaths would be reduced by 3.27% and 5.2%,
respectively; and by considering the vaccine efficacy to mutant strain as 90%,
the cumulative cases and deaths would be reduced by 0.9% and 2.5%,
respectively. Our study demonstrates that the COVID-19 pandemic may be worse
due to the occurrence of oscillations for certain mutation rates (i.e.,
outbreaks will occur repeatedly) but better due to stability at a lower
infection level with a larger mutation rate. We perform sensitivity analysis
using the Latin Hypercube Sampling (LHS) methodology and partial rank
correlation coefficients (PRCCs) to illustrate the impact of parameters on the
basic reproduction number, the number of cumulative cases, and the number of
deaths, which ultimately sheds light on disease mitigation.

Keywords: COVID-19; Endemic bubble; Mutation; Imperfect vaccination; MCMC; Two
strain dynamics; Hopf bifurcation; Transcritical bifurcation; Sensitivity
analysis.

##  1 Introduction

Contagious diseases are one of the foremost reasons for demise worldwide. The
spread of contagious diseases dangerously affects the growth of countries and
the evolution of a population. Though modern scientific medicine has made
rapid advancements, the diseases have not been completely eradicated. Diseases
have obtained new versions due to the genetic variations of pathogens
triggered via mutations. Many pathogens are characterized by more than one
variant [65, 58]. Virus or pathogens mutations are general in contagious
diseases such as HBV [65], Influenza [58], and HIV [22]. Multi-strain of 1918
avian influenza virus with the mutation have been recognized by Iwami et al.
[35]. The existence of different variants of a pathogen is mainly due to
resist immune attacks of the host or induced by treatment with antiviral drugs
or antibodies [22]. Ultimately, they confirm the persistence of disease in a
host. Sansonetti and Arondel [64] have revealed that mutant strains can be
associated with higher virulence to disease than the native strains, and those
people diseased with mutant strains have a higher death rate in the contagious
diseases such as plague, influenza A, etc. Thus, one of the major challenges
in stopping the spread of infectious diseases is to treat with the genetic
variations of pathogens [52, 59, 44]. Mathematical models are helpful to
describe and understand the dynamics of different strains under mutation.
Various epidemic models with multi-strain contacts and mutation have been
proposed in recent years from different aspects [44, 12, 43, 51]. Liu et al.
[44] proposed a mathematical model for Influenza with virus mutation and
analyzed the model in the sense of permanence of the disease. Cai et al. [12]
proposed a two-strain model with vaccination. The authors studied the
existence and stability of the equilibria as well as the existence of Hopf
bifurcation from endemic equilibria. Li et al. [43] proposed a two-strain SIR
model with infection age and mutation. The authors analyzed the stability of
equilibria and Hopf bifurcation. Moreover, epidemiological investigations have
exposed that the phenomenon of mutations leads to further resistant viruses
giving the emergence of many new dangerous epidemics or even new serious
pandemics.

Since December 2019, a novel coronavirus (identified to produce a respiratory
disease acknowledged as COVID-19) has spread broadly and rapidly and has since
oppressed a significant proportion of the population worldwide. The World
Health Organization (WHO) has detected the spread of COVID-19 as a pandemic,
and as of December 4, 2021, over 263 million peoples were diseased, and about
5.2 million died of the virus. The SARS-CoV-2 virus triggered by severe acute
respiratory syndrome is also mutating. Lately, numerous variants of the SARS-
CoV-2 virus have been identified. These variants are described according to
the number and types of mutations [39, 42]. Since the start of the COVID-19
pandemic, the SARS-CoV-2 coronavirus that initiated COVID-19 has mutated,
rising in different variants of the virus. Numerous SARS-CoV-2 variants have
developed worldwide, and the presence of different variants depends on several
factors. One of these is called the delta variant, which was first identified
in India [17]. Different variants have appeared in Brazil, England,
California, and other countries. More transmittable variants such as beta,
which first emerged in South Africa, may have improved the ability to re-
infect individuals who have recuperated from previous versions of the virus
and also be somewhat resistant to some of the coronavirus vaccines in
development [36]. These new variants might have distinct features that can
influence the death rate and transmissibility [39, 42, 25]. From October 2020,
the number of infected cases of SARS-CoV-2 and related deaths augmented
drastically in England. It has been discovered that the new SARS-CoV-2 variant
VOC-202012/01 was widespread, and its proportion amplified throughout the
latest months in England [10, 62]. The mutations of viruses are frequent, and,
as an outcome, SARS-CoV-2 can develop mutations with immunological resistance
and fitness advantages [39]. It is anticipated that further mutations will
occur worldwide and probably even more after worldwide vaccination due to
mutation force [63]. Therefore, analyzing the effect of new strains of the
SARS-CoV-2 virus is supremely significant. In the literature, a few authors
proposed multi-variant mathematical models for COVID-19 [25, 38, 6]. Gonzalez-
Parra et al. [25] studied the effect of a new, more infectious SARS-CoV-2
variant (VOC-202012/01 of lineage B.1.1.7) on hospitalizations, prevalence,
and deaths associated with the SARS-CoV-2 virus. Khyar et al. [38] proposed a
multi-strain SEIR model with general incident rate and studied the global
dynamics of the model. The authors also discussed the quarantine strategy for
controlling the disease spread and fit the model to the Moroccan clinical data
of COVID-19. Arruda et al. [6] proposed a model for COVID-19 and studied the
optimal control of multi-strain epidemics.

Vaccination has been an effective strategy in battling the spread of
contagious diseases, e.g., measles, influenza, and pertussis. In history, the
elimination of smallpox has been counted as the most notable victory of
vaccination ever recorded [72]. Several authors in multiple papers have
investigated the role of vaccination [12, 51, 14, 66, 47, 5, 27, 13, 2]. Under
vaccination concern, certain mutant strains will finally have the competitive
benefits amongst their contacts [66]. The impact of vaccination on the growth
of strain contacts in multi-strain viruses has also been analyzed in many
papers [28, 53, 54, 61, 60]. With the latest development of anti-COVID
vaccines, numerous models have been proposed to provide insight into the
impact of vaccination of a certain fraction of the populace on the dynamics of
the COVID-19 pandemic. For instance, Fudolig et al. [24] proposed a multi-
strain model with vaccination for COVID-19 and studied the local stability of
equilibria. Furthermore, there is some evidence of COVID-19’s vaccine efficacy
in Australia that mRNA vaccine has over 90 % efficacy against COVID-19
infection, ChAdOx1 nCoV-19 has efficacy of 62% against symptomatic infection
in the intended two-dose schedule, BNT162b2 mRNA vaccine, which has 95%
efficacy against symptomatic infection [48]. Also, the Pfizer-BioNTech
BNT162b2 mRNA vaccine has an efficacy over 95%, Johnson & Johnson [J&J] Ad26
has an efficacy over 67%, the AstraZeneca–Oxford ChAdOx1 nCov-19 vaccine has
an efficacy over 67%, and the Gamaleya GamCovidVac [Sputnik V] vaccine has the
efficacy over 90% [57]. According to the Ministry of Health and Family
Welfare, Government of India, the Indian vaccine also has vaccine efficacy
over 70-90% [55].

To determine effectual countermeasures, it is significant to develop
mathematical models that support us in predicting and understanding the spread
of COVID-19 and providing strategies on what could be implemented to limit its
spread. Mathematical modeling in epidemiology provides a progressively greater
room to public health research. This research discipline participates to
sufficiently comprehend the studied epidemiological phenomenon and capture the
distinct issues that can give rise to a terrible epidemic or even an alarming
pandemic worldwide. The classical susceptible-infected-recovered (SIR)
epidemic model was first proposed by Kermack et al. [37]. To obtain a sharper
understanding of various vaccination strategies and their impacts on the
number of infected individuals, Kermack-McKendrick type models have grabbed a
vital role. This type of model has been benefited to comprehend vaccination
dynamics on various diseases [2]. It is essential to emphasize that nowadays,
Kermack-McKendrick kind mathematical models have helped explain COVID-19
epidemics properties worldwide. These models have been utilized to estimate
the basic reproductive number associated with the disease and various
parameters engaged in its spread. Additional use of this type of models has
been focused on proposing and assessing the impact of different control
measures categorized as NPIs. For example, the authors in papers [8, 9, 7]
proposed mathematical models for COVID-19 and analyzed the impact of NPIs on
the disease dynamics. It is of paramount significance to develop mathematical
models that can perfectly forecast the spread of COVID-19 so that the disease
can be controlled and restrictions can be securely relaxed. However, the
infection incubation period may occupy a long time interval in some cases. An
incubated person is not yet infectious in this time interval and remains
latent. Therefore, another compartment of exposed individuals should be
included in SIR, and the new model will have SEIR abbreviation [32].

Inspired by the evidence mentioned above about imperfect vaccine and mutation
of the virus, in this paper we utilize an SEIR-type mathematical model to
comprehend the dynamics of disease spread on the human populace under
imperfect vaccination and two variants of the virus. The general methodology
and mathematical model can be inferred to enhance the number of parameters and
differential equations. We incorporate the vaccination compartment to the two-
strain model to examine the effectiveness of the anti-COVID-19 vaccination,
which is currently being employed in many countries to help battle the intense
pandemic situation. It is supposed that the spread of a virus can mutate in
the host to make a second, co-circulating, mutant strain. After some period of
infection, the original strain, referred to as native strain, is converted to
a mutant strain, such that a proportion of the people infected by the original
strain are also carrying mutant strain, thus we consider mutation in our
proposed epidemic model. We intend to study the dynamical behavior of the
strains’ contacts under the vaccination scheme and investigate the impact of
parameters (vaccination proportion, mutation rate, etc.) to demonstrate how
they influence disease transmission. We also assume that recovered individuals
from native strain have 100% immunity against native strain but can be
infected by mutant strain. The objectives of the present study are the
following: (i) Construction of an epidemic model that will describe the
dynamics of mutant strain under imperfect vaccination. (ii) Investigation of
the impact of an imperfect vaccine on the dynamics of the model. (iii)
Investigation of different bifurcations with respect to variation of various
parameters. (iv) Application of the proposed model to the data of COVID-19 in
India. (v) Observation of the COVID-19 dynamics with respect to the key
parameter related to the mutation.

The remaining paper is organized as follows. Section 2 describes the proposed
model with imperfect vaccine and mutation. Section 3 discusses the non-
negativity and boundedness of the solutions of the proposed system. Section 4
represents the dynamical analysis of the proposed system including the basic
reproduction number, existence, and stability of possible equilibria,
transcritical bifurcation, and Hopf bifurcation. Section 5 describes the
implications of disease control, and we obtained explicit expression of
vaccine-induced herd immunity. Numerical evaluations have been presented in
Section 6. Cumulative cases and cumulative mortality data for COVID-19
pandemic in India have been used to parametrize the model, and the impact of
different parameters on the cumulative cases and deaths have been shown in
Section 7. Sensitivity analysis of the parameters with respect to the basic
reproduction number, cumulative cases, and cumulative deaths has been
performed in Section 8. The paper ends with a thorough discussion in Section
9.

##  2 Model formulation

We introduce a homogeneous two-strain model with imperfect vaccination. The
system starts with a population exposed to both the native (original) and
mutant (variant of original) strains of the virus. Some previous studies [25,
19] revealed that the variant is more transmissible and severe than the
original strain, and antibody neutralization is reduced in COVID-19 patients
and vaccine recipients in various countries, including the US. Nature news has
reported that the mutant strain is spreading quickly in India and has become
the dominant strain [56]. A mutation is accounted for in epidemic models
through a term that transfers individuals infected with one of the strains
into individuals infected with the other [44, 12, 50]. We assume that
vaccination is applied only to healthy individuals, so only susceptible
individuals get vaccinated. Further, we also assume that the vaccine is
imperfect, that is, the vaccinated individuals can become infected with both
native and mutant strains of the virus. That happens at reduced transmission
rates $\delta_{1}\beta_{1}$ and $\delta_{2}\beta_{2},$ where
$0\leq\delta_{1}\leq 1,0\leq\delta_{2}\leq 1$ are the reduction coefficients
of native and mutant strains, respectively. If $\delta_{1}=0,$ then vaccinated
individuals will not get infected with native strain, i.e., the vaccine is
perfect for native strain, and $\delta_{1}=1$ means vaccinated individuals get
infected just like susceptible individuals, i.e., vaccine plays no protective
role to native strain. Here $1-\delta_{1}$ describes the vaccine efficacy to
native strain. A similar scenario applies to the mutant strain, and
$1-\delta_{2}$ describes the vaccine efficacy to the mutant strain. Some
vaccinated individuals can go back to susceptible individuals due to partial
immunity. There are some shreds of evidence that there is a chance of a second
COVID-19 infection after being diagnosed with first, from which recovered [31,
67, 18]. Therefore, we assume that the recovered individuals of native strain
can also become infected via mutant strain of the virus at a reduced rate
$\delta_{3}\beta_{2}$.

The model is composed of a system of differential equations that has eight
compartments: susceptible compartment $S-$ individuals in this compartments
are healthy but can be infected by both the native and mutant strains of the
virus; vaccinated compartment $V-$ individuals that applied to vaccination,
these individuals can also become infected by both the native and mutant
strains of the virus but at lower rates; exposed compartment $E_{1}-$
individuals that are exposed to a native strain of virus; exposed compartment
$E_{2}-$ individuals that are exposed to a mutant strain of virus; infected
compartment $I_{1}-$ individuals that are infected to a native strain of
virus; infected compartment $I_{2}-$ Individuals that are infected to a mutant
strain of virus; recovered compartment $R_{1}-$ individuals that were infected
to native strain and are now immune to the native strain but not immune to the
mutant strain; recovered compartment $R_{2}-$ individuals infected with mutant
strain are now immune to both native and mutant strains and do not interact
with the remaining compartments. The biological interpretations of the
parameters involved in the model are given in Table 1. The schematic diagram
of the model is given in Figure 1. The total population size is
$N(t)=S(t)+V(t)+E_{1}(t)+E_{2}(t)+I_{1}(t)+I_{2}(t)+R_{1}(t)+R_{2}(t).$

Table 1: Biological interpretations of parameters. Parameters |  Biological interpretations
---|---
$\Lambda$ |  The recruitment rate at which new individuals enter in the susceptible population
$\beta_{1}$ |  Infection rate of the native strain
$\beta_{2}$ |  Infection rate of the mutant strain
$p$ |  Per capita vaccination rate of susceptible individuals
$1/\mu$ |  Average life expectancy of the individuals of all compartments
$1-\delta_{1}$ |  Efficacy of vaccine to native strain
$1-\delta_{2}$ |  Efficacy of vaccine to mutant strain
$\gamma$ |  Per capita rate of lost of immunity of vaccinated individuals
$a_{1}$ |  Per capita rate at which the exposed individuals of native strain become infectious
$a_{2}$ |  Per capita rate at which the exposed individuals of mutant strain become infectious
$\alpha_{1}$ |  Per capita recovery rate of native strain
$\alpha_{2}$ |  Per capita recovery rate of mutant strain
$d_{1}$ |  Per capita death rate due to native strain
$d_{2}$ |  Per capita death rate due to mutant strain
$\delta_{3}$ |  Reduction coefficient of infection after recovery
$\nu_{1}$ |  Per capita mutation rate of native strain
Figure 1: Schematic diagram of system (1).

Based on the parameters given in Table 1 and schematic diagram 1, the dynamics
of the disease transmission can be governed by the following system of
ordinary differential equations:

$\displaystyle\frac{dS}{dt}$
$\displaystyle=\Lambda-\beta_{1}I_{1}S-\beta_{2}I_{2}S-(\mu+p)S+\gamma V,$ (1)
$\displaystyle\frac{dV}{dt}$
$\displaystyle=pS-\delta_{1}\beta_{1}I_{1}V-\delta_{2}\beta_{2}I_{2}V-(\mu+\gamma)V,$
$\displaystyle\frac{dE_{1}}{dt}$
$\displaystyle=\beta_{1}(S+\delta_{1}V)I_{1}-(a_{1}+\mu)E_{1},$
$\displaystyle\frac{dE_{2}}{dt}$
$\displaystyle=\beta_{2}(S+\delta_{2}V+\delta_{3}R_{1})I_{2}-(a_{2}+\mu)E_{2},$
$\displaystyle\frac{dI_{1}}{dt}$
$\displaystyle=a_{1}E_{1}-(\alpha_{1}+\mu+d_{1}+\nu_{1})I_{1},$
$\displaystyle\frac{dI_{2}}{dt}$
$\displaystyle=a_{2}E_{2}-(\alpha_{2}+\mu+d_{2})I_{2}+\nu_{1}I_{1},$
$\displaystyle\frac{dR_{1}}{dt}$
$\displaystyle=\alpha_{1}I_{1}-\delta_{3}\beta_{2}I_{2}R_{1}-\mu R_{1},$
$\displaystyle\frac{dR_{2}}{dt}$ $\displaystyle=\alpha_{2}I_{2}-\mu R_{2},$

with the initial conditions: $S(0)>0,V(0)\geq 0,E_{1}(0)\geq 0,E_{2}(0)\geq
0,I_{1}(0)\geq 0,I_{2}(0)\geq 0,R_{1}(0)\geq 0,R_{2}(0)\geq 0.$

##  3 Well-posedness

This section summarizes the positivity and boundedness of solutions of system
(1). Positivity is significant for biologically feasible solutions of the
system while boundedness infers that solutions are finite. System (1) is given
by the following bounded planes:

$\displaystyle\frac{dS}{dt}\Bigg{|}_{S=0,V\neq 0,E_{1}\neq 0,E_{2}\neq
0,I_{1}\neq 0,I_{2}\neq 0,R_{1}\neq 0,R_{2}\neq 0}$
$\displaystyle=\Lambda+\gamma V>0,$ $\displaystyle\frac{dV}{dt}\Bigg{|}_{S\neq
0,V=0,E_{1}\neq 0,E_{2}\neq 0,I_{1}\neq 0,I_{2}\neq 0,R_{1}\neq 0,R_{2}\neq
0}$ $\displaystyle=pS\geq 0,$ $\displaystyle\frac{dE_{1}}{dt}\Bigg{|}_{S\neq
0,V\neq 0,E_{1}=0,E_{2}\neq 0,I_{1}\neq 0,I_{2}\neq 0,R_{1}\neq 0,R_{2}\neq
0}$ $\displaystyle=\beta_{1}(S+\delta_{1}V)I_{1}\geq 0,$
$\displaystyle\frac{dE_{2}}{dt}\Bigg{|}_{S\neq 0,V\neq 0,E_{1}\neq
0,E_{2}=0,I_{1}\neq 0,I_{2}\neq 0,R_{1}\neq 0,R_{2}\neq 0}$
$\displaystyle=\beta_{2}(S+\delta_{2}V)I_{2}+\delta_{3}\beta_{2}I_{2}R_{1}\geq
0,$ $\displaystyle\frac{dI_{1}}{dt}\Bigg{|}_{S\neq 0,V\neq 0,E_{1}\neq
0,E_{2}\neq 0,I_{1}=0,I_{2}\neq 0,R_{1}\neq 0,R_{2}\neq 0}$
$\displaystyle=a_{1}E_{1}\geq 0,$
$\displaystyle\frac{dI_{2}}{dt}\Bigg{|}_{S\neq 0,V\neq 0,E_{1}\neq 0,E_{2}\neq
0,I_{1}\neq 0,I_{2}=0,R_{1}\neq 0,R_{2}\neq 0}$
$\displaystyle=a_{2}E_{2}+\nu_{1}I_{1}\geq 0,$
$\displaystyle\frac{dR_{1}}{dt}\Bigg{|}_{S\neq 0,V\neq 0,E_{1}\neq 0,E_{2}\neq
0,I_{1}\neq 0,I_{2}\neq 0,R_{1}=0,R_{2}\neq 0}$
$\displaystyle=\alpha_{1}I_{1}\geq 0,$
$\displaystyle\frac{dR_{2}}{dt}\Bigg{|}_{S\neq 0,V\neq 0,E_{1}\neq 0,E_{2}\neq
0,I_{1}\neq 0,I_{2}\neq 0,R_{1}\neq 0,R_{2}=0}$
$\displaystyle=\alpha_{2}I_{2}\geq 0.$

Note that on each of the bounding planes of the non-negative cone of
$\mathbb{R}^{8}_{+}$, all rates in the system (1) are non-negative. Thus, if
we initiate this in the interior of this cone, we shall always remain in this
cone as the direction of the vector field is inward on all the bounding
planes. Therefore, the non-negativity of all solutions is guaranteed if we
start from a non-negative initial point. Furthermore, system (1) also states
that the population $N$ follows the below differential equation:

$\frac{dN}{dt}=\Lambda-\mu N-d_{1}I_{1}-d_{2}I_{2},$ (2)

which gives

$\Lambda-(\mu+d_{1}+d_{2})N\leq\frac{dN}{dt}\leq\Lambda-\mu N.$

Now integrating the above inequality and using initial conditions, we obtain

$\frac{\Lambda}{\mu+d_{1}+d_{2}}+\left(N(0)-\frac{\Lambda}{\mu+d_{1}+d_{2}}\right)e^{-(\mu+d_{1}+d_{2})t}\leq
N(t)\leq\frac{\Lambda}{\mu}+\left(N(0)-\frac{\Lambda}{\mu}\right)e^{-\mu t},$

Considering $t\rightarrow+\infty,$ we obtain

$\frac{\Lambda}{\mu+d_{1}+d_{2}}\leq\liminf\limits_{t\rightarrow+\infty}N(t)\leq\limsup\limits_{t\rightarrow+\infty}N(t)\leq\frac{\Lambda}{\mu}.$

Hence, the feasible region for the system (1) is

$\Delta=\left\\{(S,V,E_{1},E_{2},I_{1},I_{2},R_{1},R_{2})\in\mathbb{R}^{8}_{+}:0<S+V+E_{1}+E_{2}+I_{1}+I_{2}+R_{1}+R_{2}\leq\frac{\Lambda}{\mu}\subset\mathbb{R}^{8}_{+}\right\\}.$

From the above analysis, we conclude the following consequence:

######  Theorem 3.1.

The region $\Delta$ is positively invariant and attracting.

Therefore, the system (1) is well-posed and epidemiologically feasible since
all variables remain nonnegative for all $t\geq 0$. Further, since the
equations of system (1) are continuous and have continuous partial
derivatives, then they satisfy the Lipschitz condition. Additionally, from
Theorem 3.1, system (1) is uniformly bounded. Hence, the solution of the
system (1) exists and is unique.

##  4 Rigorous analysis

This section is devoted to investigating the dynamical behaviors of system (1)
including the computation of the basic reproduction number, the existence of
possible equilibria and their stability, and possible bifurcations.

###  4.1 Disease free equilibrium (DFE) and basic reproduction number

The disease free equilibrium (DFE) can be obtained by setting all infected
variables ($E_{1},E_{2},I_{1},I_{2}$) equal to zero while all non-infected
variables ($S,V,R_{1},R_{2}$) are non-zero. The DFE of system (1) is given by

$D^{0}=(S^{0},V^{0},E_{1}^{0},E_{2}^{0},I_{1}^{0},I_{2}^{0},R_{1}^{0},R_{2}^{0})=\left(\frac{\Lambda(\mu+\gamma)}{\mu(\mu+\gamma+p)},\frac{\Lambda
p}{\mu(\mu+\gamma+p)},0,0,0,0,0,0\right).$ (3)

To obtain the basic reproduction number, we use the next generation method
[20, 68]. By considering $x=(E_{1},E_{2},I_{1},I_{2})^{T},$ we have

$x^{\prime}=f(x)-v(x),$ (4)

where

$\displaystyle f=\left(\begin{array}[]{c}\beta_{1}(S+\delta_{1}V)I_{1}\\\
\beta_{2}(S+\delta_{2}V)I_{2}+\delta_{3}\beta_{2}I_{2}R_{1}\\\ 0\\\
0\end{array}\right),\quad v=\left(\begin{array}[]{c}(a_{1}+\mu)E_{1}\\\
(a_{2}+\mu)E_{2}\\\ (\alpha_{1}+\mu+d_{1}+\nu_{1})I_{1}-a_{1}E_{1}\\\
(\alpha_{2}+\mu+d_{2})I_{2}-\nu_{1}I_{1}-a_{2}E_{2}\end{array}\right).$

The Jacobian of $f(x)$ and $v(x)$ at $D^{0}$ are

$\displaystyle
F=Df(D^{0})=\left(\begin{array}[]{cccc}0&0&\beta_{1}(S^{0}+\delta_{1}V^{0})&0\\\
0&0&0&\beta_{2}(S^{0}+\delta_{2}V^{0})\\\ 0&0&0&0\\\
0&0&0&0\end{array}\right),$

and

$\displaystyle V$
$\displaystyle=Dv(D^{0})=\left(\begin{array}[]{cccc}(a_{1}+\mu)&0&0&0\\\
0&(a_{2}+\mu)&0&0\\\ -a_{1}&0&(\alpha_{1}+\mu+d_{1}+\nu_{1})&0\\\
0&-a_{2}&-\nu_{1}&(\alpha_{2}+\mu+d_{2})\end{array}\right),$

respectively. Hence, the basic reproduction number for system (1) is the
dominant eigenvalue or spectral radius of the next generation matrix $FV^{-1}$
which is given by the following expression:

$\displaystyle R_{0}=\rho(FV^{-1})=\max\left\\{R_{1},R_{2}\right\\},$ (5)

where

$\displaystyle R_{1}$
$\displaystyle=\frac{\beta_{1}a_{1}(S^{0}+\delta_{1}V^{0})}{(\mu+a_{1})(\mu+d_{1}+\alpha_{1}+\nu_{1})}$
$\displaystyle=\frac{\Lambda\beta_{1}a_{1}(\gamma+\mu+p\delta_{1})}{\mu(\mu+a_{1})(\mu+d_{1}+\alpha_{1}+\nu_{1})(p+\gamma+\mu)},$
$\displaystyle R_{2}$
$\displaystyle=\frac{\beta_{2}a_{2}(S^{0}+\delta_{2}V^{0})}{(\mu+a_{2})(\mu+d_{2}+\alpha_{2})}$
$\displaystyle=\frac{\Lambda\beta_{2}a_{2}(\gamma+\mu+p\delta_{2})}{\mu(\mu+a_{2})(\mu+d_{2}+\alpha_{2})(p+\gamma+\mu)}.$

Here $R_{1}(R_{2})$ represents the average number of secondary infection cases
generated by a single infectious individuals of the native (mutant) strain of
the virus, called the basic reproduction number of the native (mutant) strain.

####  4.1.1 Interpretation of the basic reproduction number

As stated above, the basic reproduction number $R_{0}$ is the maximum of the
two basic reproduction numbers, $R_{1}$ and $R_{2}$. The basic reproduction
number $R_{1}$ is given by the product of the infection rate of the
susceptible (unvaccinated) and vaccinated individuals by native strain
infectious individuals (near the disease-free equilibrium),
[$\beta_{1}(S^{0}+\delta_{1}V^{0})$], the proportion of the exposed
individuals to the native strain that survived in the exposed class ($E_{1}$)
and moved to infected compartment ($I_{1}$) [$\frac{a_{1}}{\mu+a_{1}}$], and
the average time duration in the infectious class ($I_{1}$)
[$\frac{1}{\mu+d_{1}+\alpha_{1}+\nu_{1}}$]. The basic reproduction number
$R_{2}$ is given by the product of the infection rate of the susceptible
(unvaccinated) and vaccinated individuals by the mutant strain infectious
individuals (near the disease-free equilibrium),
[$\beta_{2}(S^{0}+\delta_{1}V^{0})$], the proportion of the exposed
individuals to the mutant strain that survived in the exposed class ($E_{2}$)
and moved to infected compartment ($I_{2}$) [$\frac{a_{2}}{\mu+a_{2}}$], and
the average time duration in the infectious class ($I_{2}$)
[$\frac{1}{\mu+d_{2}+\alpha_{2}}$].

######  Remark 4.1.

If $\delta_{1}=0,\delta_{2}=0,$ then vaccine is perfectly effective to both
native and mutant strains and if $p=0$ then the system reduces without
vaccination. For this case, the basic reproduction numbers of the native
strain and mutant strain, $R_{1wv}$ and $R_{2wv}$ are given by, respectively,

$\displaystyle R_{1wv}$
$\displaystyle=\frac{\Lambda\beta_{1}a_{1}}{\mu(\mu+a_{1})(\mu+d_{1}+\alpha_{1}+\nu_{1})},$
$\displaystyle R_{2wv}$
$\displaystyle=\frac{\Lambda\beta_{2}a_{2}}{\mu(\mu+a_{2})(\mu+d_{2}+\alpha_{2})}.$

Thus, the basic reproduction number of system without vaccination is given by
$R_{0wv}=\max\left\\{R_{1wv},R_{2wv}\right\\}.$

######  Theorem 4.2.

If $R_{0}=\max\left\\{R_{1},R_{2}\right\\}<1,$ the DFE $(D^{0})$ is locally
asymptotically stable; if $R_{0}=\max\left\\{R_{1},R_{2}\right\\}>1,$ $D^{0}$
is unstable.

###### Proof.

By linearizing the system (1) at $D^{0}$, we obtain the following
characteristic equation

$\displaystyle(\lambda+\mu)^{2}(\lambda+p+\gamma+\mu)(\lambda^{2}+\lambda(2\mu+a_{1}+d_{1}+\alpha_{1})+(\mu+a_{1})(\mu+d_{1}+\alpha_{1}+\nu_{1})(1-R_{1}))$
(6)
$\displaystyle(\lambda^{2}+\lambda(2\mu+a_{2}+d_{2}+\alpha_{2})+(\mu+a_{2})(\mu+d_{2}+\alpha_{2})(1-R_{2}))=0.$

From the above characteristic equation, it is easy to see that all the roots
of Eq. (6) are negative or have negative real parts for $R_{0}<1.$ Hence, the
$D^{0}$ of system (1) is locally asymptotically stable for $R_{0}<1.$ If
$R_{0}>1,$ at least one of the roots of Eq. (6) has positive real part. Hence,
$D^{0}$ is unstable for $R_{0}>1.$ ∎

######  Theorem 4.3.

$D^{0}$ is globally asymptotically stable, whenever
$R_{0}=\max\left\\{R_{1},R_{2}\right\\}<1$.

###### Proof.

To prove the global stability of $D^{0}$, we follow the approach given by
Castillo-Chavez et al. [15]. We rewrite the system (1) as follows

$\displaystyle\frac{dX}{dt}$ $\displaystyle=F(X,Y),$ (7)
$\displaystyle\frac{dY}{dt}$ $\displaystyle=G(X,Y),\quad G(X,0)=0,$

where $X=(S,V,R_{1},R_{2})\in\mathbb{R}^{4}$ signifies the number of
uninfected individuals and $Y=(E_{1},E_{2},I_{1},I_{2})\in\mathbb{R}^{4}$
signifies the number of infected individuals. Disease-free equilibrium
$(D^{0})$ is globally stable if the following two conditions are fulfilled:

  * (H1)

For $\frac{dX}{dt}=F(X,Y),$ $X^{*}$ is globally asymptotically stable,

  * (H2)

$G(X,Y)=MY-\hat{G}(X,Y),\leavevmode\nobreak\ \hat{G}(X,Y)>0$ for
$(X,Y)\in\Delta$,

where $M=D_{Y}G(X^{*},0)$ is an $M$-matrix. For the system (1), we have

$F(X,0)=\left(\begin{array}[]{c}\Lambda-(\mu+p)S+\gamma V\\\
pS-(\mu+\gamma)V\\\ 0\\\ 0\end{array}\right).$ (8)

It is obvious that the equilibrium
$X^{*}=\left(\frac{\Lambda(\mu+\gamma)}{\mu(\mu+\gamma+p)},\frac{\Lambda
p}{\mu(\mu+\gamma+p)},0,0\right)$ is globally asymptotically stable of system
(8). Further, for system (1), we obtain

$\displaystyle M$
$\displaystyle=\left(\begin{array}[]{cccc}-(\mu+a_{1})&0&\beta_{1}(S^{0}+\delta_{1}V^{0})&0\\\
0&-(\mu+a_{2})&0&\beta_{2}(S^{0}+\delta_{2}V^{0})\\\
a_{1}&0&-(\alpha_{1}+\mu+d_{1}+\nu_{1})&0\\\
0&a_{2}&\nu_{1}&-(\alpha_{2}+\mu+d_{2})\end{array}\right),$
$\displaystyle\hat{G}(X,Y)$
$\displaystyle=\left(\begin{array}[]{c}\beta_{1}I_{1}((S^{0}+\delta_{1}V^{0})-(S+\delta_{1}V))\\\
\beta_{2}I_{2}((S^{0}+\delta_{2}V^{0})-(S+\delta_{2}V+\delta_{3}R_{1}))\\\
0\\\ 0\end{array}\right).$

It is clear that $\hat{G}(X,Y)\geq 0$. Hence, $D^{0}$ is globally stable, i.e.
every solution of the system (1) approaches the DFE ($D^{0}$) as
$t\rightarrow\infty$ for $R_{0}<1.$ Thus, the disease (i.e., both native and
mutant strains) will be eliminated from the community if $R_{0}<1.$
Consequently, $R_{0}$ represents the threshold value for the existence of
other positive equilibria of the system (1). Moreover, it is eminent that the
basic reproduction number $(R_{0})$ represents the average number of secondary
infections that occurred from a single infected individual in the whole
susceptible population. Therefore, if $R_{0}<1$, each infected individual in
the entire infectious period will produce less than one infected individual on
average, which implies that the disease will die out. However, if $R_{0}>1,$
then each infected individual in the whole infectious period having contact
with susceptible individuals will generate more than one infected individual;
this leads to the disease invading the susceptible population. ∎

It must be mentioned that, for mathematical (endemic) models such as (1), the
epidemiological necessity $R_{0}<1$ is sufficient as well as necessary for
eradication of the disease. This is because, for such mathematical models
(i.e., Kermack-McKendrick models with demographic dynamics), the disease will
persist whenever $R_{0}>1$ (this is because the pool of new susceptible
individuals will continuously be refilled, by immigration or birth, thereby
letting the disease to maintain itself in the community). If the demographic
effects are not allowed (i.e., in the case of a single outbreak/epidemic model
is used), the epidemiological condition $R_{0}<1$ is only sufficient, but not
necessary, for eradicating the epidemic. For such epidemic models (with no
demographic dynamics), the disease always dies out with time (irrespective of
the value of the basic reproduction number of the epidemic models). In other
words, even if the basic reproduction number exceeds unity, the disease will
eventually die out; this is because the endemic rises and reach a peak.

###  4.2 Mutant dominant equilibrium and its stability

First of all, it should be mentioned that in the absence of native strain
($I_{1}=0$), the system (1) reduces to the following subsystem:

$\displaystyle\frac{dS}{dt}$
$\displaystyle=\Lambda-\beta_{2}I_{2}S-(\mu+p)S+\gamma V,$ (9)
$\displaystyle\frac{dV}{dt}$
$\displaystyle=pS-\delta_{2}\beta_{2}I_{2}V-(\mu+\gamma)V,$
$\displaystyle\frac{dE_{2}}{dt}$
$\displaystyle=\beta_{2}(S+\delta_{2}V)I_{2}-(a_{2}+\mu)E_{2},$
$\displaystyle\frac{dI_{2}}{dt}$
$\displaystyle=a_{2}E_{2}-(\alpha_{2}+\mu+d_{2})I_{2},$
$\displaystyle\frac{dR_{2}}{dt}$ $\displaystyle=\alpha_{2}I_{2}-\mu R_{2}.$

The analysis of the sub-system (9) will be considered in the following
positively invariant region

$\Delta_{I_{2}}=\left\\{(S,V,E_{2},I_{2},R_{2})\in\mathbb{R}^{5}_{+}:0<S+V+E_{2}+I_{2}+R_{2}\leq\frac{\Lambda}{\mu}\subset\mathbb{R}^{5}_{+}\right\\}.$

The mutant dominant equilibrium is given by
$D^{2}=(S^{2},V^{2},0,E_{2}^{2},0,I_{2}^{2},0,R_{2}^{2})$, where the
components of $D^{2}$ can be obtained by solving the equations of right hand
side of the subsystem (9). Thus,

$\displaystyle S^{2}$
$\displaystyle=\frac{\Lambda}{\mu}-\frac{p\left\\{\Lambda
a_{2}-(a_{2}+\mu)(\alpha_{2}+\mu+d_{2})I_{2}^{2}\right\\}}{\mu
a_{2}(\delta_{2}\beta_{2}I_{2}^{2}+\mu+\gamma+p)}-\frac{(a_{2}+\mu)(\alpha_{2}+\mu+d_{2})I_{2}^{2}}{\mu
a_{2}},$ $\displaystyle V^{2}$ $\displaystyle=\frac{p\left\\{\Lambda
a_{2}-(a_{2}+\mu)(\alpha_{2}+\mu+d_{2})I_{2}^{2}\right\\}}{\mu
a_{2}(\delta_{2}\beta_{2}I_{2}^{2}+\mu+\gamma+p)},\quad
E_{2}^{2}=\frac{(\alpha_{2}+\mu+d_{2})I_{2}^{2}}{a_{2}},\quad
R_{2}^{2}=\frac{\alpha_{2}I_{2}^{2}}{\mu},$

and $I_{2}^{2}$ satisfies the following equation:

$\displaystyle
k_{1}^{{}^{\prime}}I_{2}^{2}+k_{2}^{{}^{\prime}}I_{2}+k_{3}^{{}^{\prime}}=0,$
(10)

where

$\displaystyle k_{1}^{{}^{\prime}}$
$\displaystyle=\beta_{2}^{2}(\mu+a_{2})(\mu+d_{2}+\alpha_{2})\delta_{2},$
$\displaystyle k_{2}^{{}^{\prime}}$
$\displaystyle=\beta_{2}((\mu+a_{2})(\mu+d_{2}+\alpha_{2})(\gamma+\mu+(p+\mu)\delta_{2})-\Lambda\beta_{2}a_{2}\delta_{2}),$
$\displaystyle k_{3}^{{}^{\prime}}$
$\displaystyle=\mu(\mu+a_{2})(\mu+d_{2}+\alpha_{2})(\gamma+\mu+p)(1-R_{2}).$

We can see that Eq. (10) has zero, one, or two roots, depending on parameter
values. For the case $0<\delta_{2}\leq 1,$ $k_{3}^{{}^{\prime}}<0$ if
$R_{2}>1$, and $k_{3}^{{}^{\prime}}>0$ if $R_{2}<1$. Since Eq. (10) is
quadratic equation, therefore if $R_{2}>1$, then Eq. (10) has a unique
positive root and there is a unique mutant dominant equilibrium. If $R_{2}=1$,
then $k_{3}^{{}^{\prime}}=0$ and there is unique non-zero solution of (10),
given by $I_{2}^{2}=-\frac{k_{2}^{{}^{\prime}}}{k_{1}^{{}^{\prime}}},$ which
is positive if and only if $k_{2}^{{}^{\prime}}<0.$ If
$R_{2}=1,k_{3}^{{}^{\prime}}=0,$ then

$\Lambda\beta_{2}a_{2}(\gamma+\mu+p\delta_{2})=\mu(\mu+a_{2})(\mu+d_{2}+\alpha_{2})(p+\gamma+\mu).$
(11)

The condition $k_{2}^{{}^{\prime}}<0$ gives

$(\mu+a_{2})(\mu+d_{2}+\alpha_{2})(\gamma+\mu+(p+\mu)\delta_{2})<\Lambda\beta_{2}a_{2}\delta_{2},$

combined with (11), we obtain

$(\gamma+\mu)^{2}+(p\delta_{2})^{2}+\mu p\delta_{2}^{2}+2\gamma
p\delta_{2}+\mu p\delta_{2}<0,$

which is not possible. Hence, if $R_{2}\leq 1$, system (1) has no mutant
dominant equilibrium. Furthermore, it should be stated that for the
equilibrium $D^{2}$ to exist, it is necessary that the native strain dies out
asymptotically (i.e., $R_{1}\leq 1$). Thus, we conclude that system (9) has a
unique mutant dominant equilibrium ($D^{2}$) whenever $R_{2}>1$ and $R_{1}\leq
1.$ Further, for the stability of the equilibrium $D^{2}$, we have the
following result:

######  Theorem 4.4.

The unique mutant dominant equilibrium ($D^{2}$) is globally asymptotically
stable whenever $R_{2}>1$ and $R_{1}\leq 1.$

###### Proof.

We consider

$\displaystyle x=\frac{S}{S^{2}},\leavevmode\nobreak\
y=\frac{V}{V^{2}},\leavevmode\nobreak\
z=\frac{E_{2}}{E_{2}^{2}},\leavevmode\nobreak\ u=\frac{I_{2}}{I_{2}^{2}},$

and with the help of Eqs. of right hand side of the system (9), the system (9)
can be rewritten as follows:

$\displaystyle x^{\prime}$
$\displaystyle=x\Big{[}\frac{\Lambda}{S^{2}}\Big{(}\frac{1}{x}-1\Big{)}-\beta_{2}I_{2}^{2}(u-1)+\frac{\gamma
V^{2}}{S^{2}}\Big{(}\frac{y}{x}-1\Big{)}\Big{]},$ (12) $\displaystyle
y^{\prime}$
$\displaystyle=y\Big{[}\frac{pS^{2}}{V^{2}}\Big{(}\frac{x}{y}-1\Big{)}-\delta_{2}\beta_{2}I_{2}^{2}(u-1)\Big{]},$
$\displaystyle z^{\prime}$
$\displaystyle=z\frac{\beta_{2}I_{2}^{2}}{E_{2}^{2}}\Big{[}S^{2}\Big{(}\frac{xu}{z}-1\Big{)}+\delta_{2}V^{2}\Big{(}\frac{yu}{z}-1\Big{)}\Big{]},$
$\displaystyle u^{\prime}$
$\displaystyle=u\frac{a_{2}E_{2}^{2}}{I_{2}^{2}}\Big{[}\frac{z}{u}-1\Big{]}.$

Further, we consider the following Lyapunov function

$\displaystyle Z=k_{1}S^{2}(x-1-\ln x)+k_{2}V^{2}(y-1-\ln
y)+k_{3}E_{2}^{2}(z-1-\ln z)+k_{4}I_{2}^{2}(u-1-\ln u),$

where the positive constants $k_{1},k_{2},k_{3},$ and $k_{4}$ will be defined
below. By differentiating $Z$ with respect to $t$ along the solutions of (9),
we obtain

$\displaystyle Z^{\prime}=$ $\displaystyle
k_{1}(x-1)\Big{[}\Lambda\Big{(}\frac{1}{x}-1\Big{)}-\beta_{2}S^{2}I_{2}^{2}(u-1)+\gamma
V^{2}\Big{(}\frac{y}{x}-1\Big{)}\Big{]}$
$\displaystyle+k_{2}(y-1)\Big{[}ps^{2}\Big{(}\frac{x}{y}-1\Big{)}-\delta_{2}\beta_{2}V^{2}I_{2}^{2}(u-1)\Big{]}$
$\displaystyle+k_{3}\beta_{2}I_{2}^{2}(z-1)\Big{[}S^{2}\Big{(}\frac{xu}{z}-1\Big{)}+\delta_{2}V^{2}\Big{(}\frac{yu}{z}-1\Big{)}\Big{]}+k_{4}a_{2}E_{2}^{2}(u-1)\Big{(}\frac{z}{u}-1\Big{)}$
$\displaystyle=$ $\displaystyle k_{1}(2\Lambda+\gamma
V^{2}-\beta_{2}S^{2}I_{2}^{2})+k_{2}(pS^{2}-\delta_{2}\beta_{2}V^{2}I_{2}^{2})+k_{3}(\beta_{2}S^{2}I_{2}^{2}+\delta_{2}\beta_{2}V^{2}I_{2}^{2})$
$\displaystyle+k_{4}a_{2}E_{2}^{2}-(k_{1}\Lambda+k_{1}\gamma
V^{2}-k_{1}\beta_{2}S^{2}I_{2}^{2}-k_{2}pS^{2})x-k_{1}\Lambda\frac{1}{x}$
$\displaystyle-(-k_{1}\gamma
V^{2}+k_{2}pS^{2}-k_{2}\delta_{2}\beta_{2}V^{2}I_{2}^{2})y-k_{1}\gamma
V^{2}\frac{y}{x}-k_{2}pS^{2}\frac{x}{y}$
$\displaystyle-(k_{1}\beta_{2}S^{2}I_{2}^{2}-k_{3}\beta_{2}S^{2}I_{2}^{2})ux-(-k_{1}\beta_{2}S^{2}I_{2}^{2}-k_{2}\delta_{2}\beta_{2}V^{2}I_{2}^{2}+k_{4}a_{2}E_{2}^{2})u$
$\displaystyle-(k_{2}\delta_{2}\beta_{2}V^{2}I_{2}^{2}-k_{3}\delta_{2}\beta_{2}V^{2}I_{2}^{2})yu-(k_{3}\beta_{2}S^{2}I_{2}^{2}+k_{3}\delta_{2}\beta_{2}V^{2}I_{2}^{2}-k_{4}a_{2}E_{2}^{2})z$
$\displaystyle-
k_{3}\beta_{2}S^{2}I_{2}^{2}\frac{xu}{z}-k_{3}\delta_{2}\beta_{2}V^{2}I_{2}^{2}\frac{yu}{z}-k_{4}a_{2}E_{2}^{2}\frac{z}{u}$
$\displaystyle=:$ $\displaystyle G(x,y,z,u).$

Choose the positive constants $k_{1},k_{2},k_{3},$ and $k_{4}$ as follows:

$k_{1}=k_{2}=k_{3}=1,k_{4}=\frac{a_{2}+\mu}{a_{2}}.$

Substituting the above values into the function $G(x,y,z,u)$ gives

$\displaystyle G(x,y,z,u)=$ $\displaystyle 2\Lambda+\gamma
V^{2}+pS^{2}+(a_{2}+\mu)E_{2}^{2}-\mu S^{2}x-\Lambda\frac{1}{x}-\mu
V^{2}y-\gamma V^{2}\frac{y}{x}$ $\displaystyle-
pS^{2}\frac{x}{y}-\beta_{2}S^{2}I_{2}^{2}\frac{xu}{z}-\delta_{2}\beta_{2}V^{2}I_{2}^{2}\frac{yu}{z}-(a_{2}+\mu)E_{2}^{2}\frac{z}{u}$
$\displaystyle=$ $\displaystyle\mu S^{2}\Big{(}2-x-\frac{1}{x}\Big{)}+\gamma
V^{2}\Big{(}2-\frac{x}{y}-\frac{y}{x}\Big{)}+\mu
V^{2}\Big{(}3-\frac{1}{x}-y-\frac{x}{y}\Big{)}$
$\displaystyle+\beta_{2}S^{2}I_{2}^{2}\Big{(}3-\frac{1}{x}-\frac{xu}{z}-\frac{z}{u}\Big{)}+\delta_{2}\beta_{2}V^{2}I_{2}^{2}\Big{(}4-\frac{1}{x}-\frac{x}{y}-\frac{yu}{z}-\frac{z}{u}\Big{)}.$

By the property that the arithmetic mean is greater than or equal to the
geometric mean, $G(x,y,z,u)\leq 0,$ and the equality holds only for $x=y=1$
and $z=u,$ i.e.,

$\left\\{(x,y,z,u)\in\Delta_{I_{2}}:G(x,y,z,u)=0\right\\}\equiv\left\\{(x,y,z,u):x=y=1,z=u\right\\},$

which corresponds to the set

$\Delta_{I_{2}}^{\prime}=\left\\{(S,V,E_{2},I_{2}):S=S^{*},V=V^{*},\frac{E_{2}}{E_{2}^{*}}=\frac{I_{2}}{I_{2}^{*}}\right\\}\subset\Delta_{I_{2}}\subset\Delta.$

It is evident to see that the maximum invariant set of (9) on the set
$\Delta_{I_{2}}^{\prime}$ is the singleton $\left\\{D^{2}\right\\}$, then the
mutant dominant equilibrium $D^{2}$ is globally stable in
$\Delta_{I_{2}}\subset\Delta$ by LaSalle’s Invariance Principle [40]. ∎

###  4.3 Coexistence equilibrium and its stability

This section examines the existence and global stability of the coexistence
equilibrium of system (1). First of all, let us consider the endemic
equilibrium
$D^{*}=(S^{*},V^{*},E_{1}^{*},E_{2}^{*},I_{1}^{*},I_{2}^{*},R_{1}^{*},R_{2}^{*})$,
then $S^{*},V^{*},$ $E_{1}^{*},E_{2}^{*},I_{1}^{*},I_{2}^{*},R_{1}^{*},$ and
$R_{2}^{*}$ satisfy the following equations:

$\displaystyle\Lambda-\beta_{1}I_{1}S-\beta_{2}I_{2}S-(\mu+p)S+\gamma V=0,$
(13) $\displaystyle
pS-\delta_{1}\beta_{1}I_{1}V-\delta_{2}\beta_{2}I_{2}V-(\mu+\gamma)V=0,$
$\displaystyle\beta_{1}(S+\delta_{1}V)I_{1}-(a_{1}+\mu)E_{1}=0,$
$\displaystyle\beta_{2}(S+\delta_{2}V+\delta_{3}R_{1})I_{2}-(a_{2}+\mu)E_{2}=0,$
$\displaystyle a_{1}E_{1}-(\alpha_{1}+\mu+d_{1}+\nu_{1})I_{1}=0,$
$\displaystyle a_{2}E_{2}-(\alpha_{2}+\mu+d_{2})I_{2}+\nu_{1}I_{1}=0,$
$\displaystyle\alpha_{1}I_{1}-\delta_{3}\beta_{2}I_{2}R_{1}-\mu R_{1}=0,$
$\displaystyle\alpha_{2}I_{2}-\mu R_{2}=0.$

The above Eqs (13) lead the following expressions:

$\displaystyle S^{*}$
$\displaystyle=\frac{\Lambda(\beta_{1}\delta_{1}I_{1}^{*}+\beta_{2}\delta_{2}I_{2}^{*}+\mu+\gamma)}{(\beta_{1}I_{1}^{*}+\beta_{2}I_{2}^{*}+\mu+p)(\beta_{1}\delta_{1}I_{1}^{*}+\beta_{2}\delta_{2}I_{2}^{*}+\mu+\gamma)-p\gamma},$
$\displaystyle V^{*}$ $\displaystyle=\frac{\Lambda
p}{(\beta_{1}I_{1}^{*}+\beta_{2}I_{2}^{*}+\mu+p)(\beta_{1}\delta_{1}I_{1}^{*}+\beta_{2}\delta_{2}I_{2}^{*}+\mu+\gamma)-p\gamma},$
$\displaystyle E_{1}^{*}$
$\displaystyle=\frac{(\alpha_{1}+\mu+d_{1}+\nu_{1})I_{1}^{*}}{a_{1}},\quad
E_{2}^{*}=\frac{(\alpha_{2}+\mu+d_{2})I_{2}^{*}-\nu_{1}I_{1}^{*}}{a_{2}},$
$\displaystyle R_{1}^{*}$
$\displaystyle=\frac{\alpha_{1}I_{1}^{*}}{\delta_{2}\beta_{2}I_{2}^{*}+\mu},\quad
R_{2}^{*}=\frac{\alpha_{2}I_{2}^{*}}{\mu},$

and $I_{1}^{*},I_{2}^{*}$ are the solutions of the following equations:

$\displaystyle F(I_{1},I_{2})\equiv$
$\displaystyle\beta_{2}^{2}\delta_{2}I_{2}^{2}+\beta_{2}\Big{(}-\frac{\Lambda\beta_{1}a_{1}\delta_{2}}{(\mu+a_{1})(\mu+d_{1}+\alpha_{1}+\nu_{1})}+\gamma+\mu+(p+\mu)\delta_{2}+\beta_{1}(\delta_{1}+\delta_{2})I_{1}\Big{)}I_{2}$
(14)
$\displaystyle+\beta_{1}^{2}\delta_{1}I_{1}^{2}+\beta_{1}\Big{(}-\frac{\Lambda\beta_{1}a_{1}\delta_{1}}{(\mu+a_{1})(\mu+d_{1}+\alpha_{1}+\nu_{1})}+\gamma+\mu+(p+\mu)\delta_{1}\Big{)}I_{1}$
$\displaystyle+\mu(p+\mu+\gamma)(1-R_{1})=0,$ $\displaystyle
G(I_{1},I_{2})\equiv$ $\displaystyle
I_{2}^{3}\Bigg{(}\beta_{2}^{2}\delta_{2}\Big{(}1-\frac{\Lambda
a_{2}\beta_{2}\delta_{2}}{(\mu+a_{2})(\mu+d_{2}+\alpha_{2})}-\frac{a_{2}\alpha_{1}\beta_{2}\delta_{3}I_{1}}{(\mu+a_{2})(\mu+d_{2}+\alpha_{2})}\Big{)}\Bigg{)}$
$\displaystyle+I_{2}^{2}\Bigg{(}\beta_{2}(\gamma+\mu+\delta_{2}(p+\mu))-\frac{a_{2}\beta_{2}^{2}\delta_{2}\Lambda(\gamma+2\mu-p\delta_{2})}{(\mu+a_{2})(\mu+d_{2}+\alpha_{2})}-\frac{I_{1}^{2}a_{2}\alpha_{1}\beta_{1}\beta_{2}^{2}\delta_{3}(\delta_{1}+\delta_{2})}{(\mu+a_{2})(\mu+d_{2}+\alpha_{2})}$
$\displaystyle+I_{1}\Big{(}\beta_{1}\beta_{2}(\delta_{1}+\delta_{2})-\frac{\beta_{2}^{2}(a_{2}\alpha_{1}\delta_{3}(\gamma+\mu+(p+\mu)\delta_{2})+\Lambda
a_{2}\beta_{1}\delta_{1}\delta_{2}+(\mu+a_{2})\delta_{1}\nu_{1})}{(\mu+a_{2})(\mu+d_{2}+\alpha_{2})}\Big{)}\Bigg{)}$
$\displaystyle+I_{2}\Bigg{(}\mu(p+\gamma+\mu)-\frac{\Lambda\beta_{2}a_{2}\mu(\gamma+\mu+p\delta_{2})}{(\mu+a_{2})(\mu+d_{2}+\alpha_{2})}-\frac{a_{2}\alpha_{1}\beta_{1}^{2}\beta_{2}\delta_{1}\delta_{3}I_{1}^{3}}{(\mu+a_{2})(\mu+d_{2}+\alpha_{2})}$
$\displaystyle+I_{1}^{2}\Big{(}1-\frac{\beta_{1}\beta_{2}\delta_{3}\alpha_{1}a_{2}(\gamma+\mu+(p+\mu)\delta_{1})}{(\mu+a_{2})(\mu+d_{2}+\alpha_{2})}+\frac{(\delta_{1}+\delta_{2})\nu_{1}}{(\mu+d_{2}+\alpha_{2})}\Big{)}$
$\displaystyle+I_{1}\Big{(}\frac{\nu_{1}(\gamma+\mu+(p+\mu)\delta_{2})-\beta_{2}\mu
a_{2}\alpha_{1}\delta_{3}(\gamma+\mu+p)-\Lambda\mu
a_{2}\beta_{1}\beta_{2}\delta_{1}}{(\mu+a_{2})(\mu+d_{2}+\alpha_{2})}$
$\displaystyle+\beta_{1}(\gamma+\mu+(p+\mu)\delta_{1})\Big{)}\Bigg{)}-\frac{I_{1}}{(\mu+d_{2}+\alpha_{2})}\Bigg{(}\nu_{1}\mu(p+\gamma+\mu)$
$\displaystyle+I_{1}\beta_{1}\nu_{1}(\gamma+\mu+(p+\mu)\delta_{1})+I_{1}^{2}\beta_{1}^{2}\delta_{1}\nu_{1}\Bigg{)}=0.$

If the system (14) admits a solution, then the system (1) will have an endemic
equilibrium. Obtaining the explicit expression for the exact solution of the
non-linear autonomous system (14) is a daunting task. However, we prove the
global stability of the endemic equilibrium for the special case
($\delta_{3}=0$ and $\nu_{1}=0$) of system (1) in the subsequent theorem:

######  Theorem 4.5.

If the endemic equilibrium
$D^{*}=(S^{*},V^{*},E_{1}^{*},E_{2}^{*},I_{1}^{*},I_{2}^{*},R_{1}^{*},R_{2}^{*})$
exists for $\delta_{3}=0$ and $\nu_{1}=0,$ then it is globally asymptotically
stable.

###### Proof.

For $\delta_{3}=0$ and $\nu_{1}=0,$ the endemic equilibrium
$D^{*}=(S^{*},V^{*},E_{1}^{*},E_{2}^{*},I_{1}^{*},I_{2}^{*}),$ $S^{*},V^{*},$
$E_{1}^{*},E_{2}^{*},$ $I_{1}^{*},$ and $I_{2}^{*}$ satisfy the following
equations:

$\displaystyle\Lambda-\beta_{1}I_{1}S-\beta_{2}I_{2}S-(\mu+p)S+\gamma V=0,$
(15) $\displaystyle
pS-\delta_{1}\beta_{1}I_{1}V-\delta_{2}\beta_{2}I_{2}V-(\mu+\gamma)V=0,$
$\displaystyle\beta_{1}(S+\delta_{1}V)I_{1}-(a_{1}+\mu)E_{1}=0,$
$\displaystyle\beta_{2}(S+\delta_{2}V)I_{2}-(a_{2}+\mu)E_{2}=0,$
$\displaystyle a_{1}E_{1}-(\alpha_{1}+\mu+d_{1})I_{1}=0,$ $\displaystyle
a_{2}E_{2}-(\alpha_{2}+\mu+d_{2})I_{2}=0.$

Since $R_{1}$ and $R_{2}$ do not appear explicitly in the first six equations,
therefore we omit them. By applying Eq. (15) and denoting

$x=\frac{S}{S^{*}},\leavevmode\nobreak\ y=\frac{V}{V^{*}},\leavevmode\nobreak\
z=\frac{E_{1}}{E_{1}^{*}},\leavevmode\nobreak\
u=\frac{E_{2}}{E_{2}^{*}},\leavevmode\nobreak\
v=\frac{I_{1}}{I_{1}^{*}},\leavevmode\nobreak\ w=\frac{I_{2}}{I_{2}^{*}},$

the system (1) can be rewritten as follows:

$\displaystyle x^{\prime}$
$\displaystyle=x\Big{[}\frac{\Lambda}{S^{*}}\Big{(}\frac{1}{x}-1\Big{)}-\beta_{1}I_{1}^{*}(v-1)-\beta_{2}I_{2}^{*}(w-1)+\frac{\gamma
V^{*}}{S^{*}}\Big{(}\frac{y}{x}-1\Big{)}\Big{]},$ (16) $\displaystyle
y^{\prime}$
$\displaystyle=y\Big{[}\frac{pS^{*}}{V^{*}}\Big{(}\frac{x}{y}-1\Big{)}-\delta_{1}\beta_{1}I_{1}^{*}(v-1)-\delta_{2}\beta_{2}I_{2}^{*}(w-1)\Big{]},$
$\displaystyle z^{\prime}$
$\displaystyle=z\frac{\beta_{1}I_{1}^{*}}{E_{1}^{*}}\Big{[}S^{*}\Big{(}\frac{xv}{z}-1\Big{)}+\delta_{1}V^{*}\Big{(}\frac{yv}{z}-1\Big{)}\Big{]},$
$\displaystyle u^{\prime}$
$\displaystyle=u\frac{\beta_{2}I_{2}^{*}}{E_{2}^{*}}\Big{[}S^{*}\Big{(}\frac{xw}{u}-1\Big{)}+\delta_{2}V^{*}\Big{(}\frac{yw}{u}-1\Big{)}\Big{]},$
$\displaystyle v^{\prime}$
$\displaystyle=v\frac{a_{1}E_{1}^{*}}{I_{1}^{*}}\Big{[}\frac{z}{v}-1\Big{]},$
$\displaystyle w^{\prime}$
$\displaystyle=w\frac{a_{2}E_{2}^{*}}{I_{2}^{*}}\Big{[}\frac{u}{w}-1\Big{]}.$

It is clear that the endemic equilibrium $D^{*}$ of (1) corresponds to the
positive equilibrium $\bar{D}^{*}(1,1,1,1,1,1)$ of (16), and that the global
stability of $\bar{D}^{*}$ is same as that of $D^{*},$ therefore, we will
discuss the global stability of the equilibrium $\bar{D}^{*}$ of system (16)
instead of $D^{*}.$

Define the Lyapunov function

$\displaystyle L=$ $\displaystyle k_{1}S^{*}(x-1-\ln x)+k_{2}V^{*}(y-1-\ln
y)+k_{3}E_{1}^{*}(z-1-\ln z)+k_{4}E_{2}^{*}(u-1-\ln u)$
$\displaystyle+k_{5}I_{1}^{*}(v-1-\ln v)+k_{6}I_{2}^{*}(w-1-\ln w),$

where the positive numbers $k_{1},k_{2},k_{3},k_{4},k_{5},$ and $k_{6}$ will
be given below, then differentiating $L$ with respect to $t$ along solutions
of (16), we obtain

$\displaystyle L^{\prime}=$ $\displaystyle
k_{1}(x-1)\Big{[}\Lambda\Big{(}\frac{1}{x}-1\Big{)}-\beta_{1}I_{1}^{*}S^{*}(v-1)-\beta_{2}I_{2}^{*}S^{*}(w-1)+\gamma
V^{*}\Big{(}\frac{y}{x}-1\Big{)}\Big{]}$
$\displaystyle+k_{2}(y-1)\Big{[}pS^{*}\Big{(}\frac{x}{y}-1\Big{)}-\delta_{1}\beta_{1}I_{1}^{*}V^{*}(v-1)-\delta_{2}\beta_{2}I_{2}^{*}V^{*}(w-1)\Big{]}$
$\displaystyle+k_{3}(z-1)\beta_{1}I_{1}^{*}\Big{[}S^{*}\Big{(}\frac{xv}{z}-1\Big{)}+\delta_{1}V^{*}\Big{(}\frac{yv}{z}-1\Big{)}\Big{]}$
$\displaystyle+k_{4}(u-1)\beta_{2}I_{2}^{*}\Big{[}S^{*}\Big{(}\frac{xw}{u}-1\Big{)}+\delta_{2}V^{*}\Big{(}\frac{yw}{u}-1\Big{)}\Big{]}$
$\displaystyle+k_{5}(v-1)a_{1}E_{1}^{*}\Big{[}\frac{z}{v}-1\Big{]}+k_{6}(w-1)a_{2}E_{2}^{*}\Big{[}\frac{u}{w}-1\Big{]}$
$\displaystyle=$ $\displaystyle
k_{1}(2\Lambda-\beta_{1}S^{*}I_{1}^{*}-\beta_{2}S^{*}I_{2}^{*}+\gamma
V^{*})+k_{2}(pS^{*}-\delta_{1}\beta_{1}V^{*}I_{1}^{*}-\delta_{2}\beta_{2}V^{*}I_{2}^{*})$
$\displaystyle+k_{3}(\beta_{1}S^{*}I_{1}^{*}+\delta_{1}\beta_{1}V^{*}I_{1}^{*})+k_{4}(\beta_{2}S^{*}I_{2}^{*}+\delta_{2}\beta_{2}V^{*}I_{2}^{*})+k_{5}a_{1}E_{1}^{*}+k_{6}a_{2}E_{2}^{*}$
$\displaystyle-(k_{1}\Lambda-
k_{1}\beta_{1}S^{*}I_{1}^{*}-k_{1}\beta_{2}S^{*}I_{2}^{*}+k_{1}\gamma
V^{*}-k_{2}pS^{*})x-k_{1}\Lambda\frac{1}{x}$ $\displaystyle-(-k_{1}\gamma
V^{*}+k_{2}pS^{*}-k_{2}\delta_{1}\beta_{1}V^{*}I_{1}^{*}-k_{2}\delta_{2}\beta_{2}V^{*}I_{2}^{*})y-(k_{1}\beta_{1}S^{*}I_{1}^{*}-k_{3}\beta_{1}S^{*}I_{1}^{*})xv$
$\displaystyle-(-k_{1}\beta_{1}S^{*}I_{1}^{*}-k_{2}\delta_{1}\beta_{1}V^{*}I_{1}^{*}+k_{5}a_{1}E_{1}^{*})v-(k_{1}\beta_{2}S^{*}I_{2}^{*}-k_{4}\beta_{2}S^{*}I_{2}^{*})xw$
$\displaystyle-(-k_{1}\beta_{2}S^{*}I_{2}^{*}-k_{2}\delta_{2}\beta_{2}V^{*}I_{2}^{*}+k_{6}a_{2}E_{2}^{*})w-k_{1}\gamma
V^{*}\frac{y}{x}-k_{2}pS^{*}\frac{x}{y}-k_{3}\beta_{1}S^{*}I_{1}^{*}\frac{xv}{z}$
$\displaystyle-
k_{3}\beta_{1}\delta_{1}V^{*}I_{1}^{*}\frac{yv}{z}-k_{4}\beta_{2}S^{*}I_{2}^{*}\frac{xw}{u}-k_{4}\beta_{2}\delta_{2}V^{*}I_{2}^{*}\frac{yw}{u}-(k_{1}\delta_{1}\beta_{1}V^{*}I_{1}^{*}-k_{3}\delta_{1}\beta_{1}V^{*}I_{1}^{*})yv$
$\displaystyle-(k_{2}\delta_{2}\beta_{2}V^{*}I_{2}^{*}-k_{4}\delta_{2}\beta_{2}V^{*}I_{2}^{*})yw-(k_{3}\beta_{1}S^{*}I_{1}^{*}+k_{3}\beta_{1}\delta_{1}V^{*}I_{1}^{*}-k_{5}a_{1}E_{1}^{*})z$
$\displaystyle-(k_{4}\beta_{2}S^{*}I_{2}^{*}+k_{4}\beta_{2}\delta_{2}V^{*}I_{2}^{*}-k_{6}a_{2}E_{2}^{*})u-k_{5}a_{1}E_{1}^{*}\frac{z}{v}-k_{6}a_{2}E_{2}^{*}\frac{u}{w}$
$\displaystyle=:$ $\displaystyle F(x,y,z,u,v,w).$

Now we choose the positive constants $k_{1},k_{2},k_{3},k_{4},k_{5},k_{6}$ as
follows:

$k_{1}=k_{2}=k_{3}=k_{4}=1,\leavevmode\nobreak\
k_{5}=\frac{a_{1}+\mu}{a_{1}},\leavevmode\nobreak\
k_{6}=\frac{a_{2}+\mu}{a_{2}}.$

Substituting them into the function $F(x,y,z,u,v,w)$ gives

$\displaystyle F(x,y,z,u,v,w)=$ $\displaystyle[2\Lambda+\gamma
V^{*}+pS^{*}+(a_{1}+\mu)E_{1}^{*}+(a_{2}+\mu)E_{2}^{*}]-\mu
S^{*}x-\Lambda\frac{1}{x}-\mu V^{*}y$ $\displaystyle-\gamma
V^{*}\frac{y}{x}-pS^{*}\frac{x}{y}-\beta_{1}S^{*}I_{1}^{*}\frac{xv}{z}-\delta_{1}\beta_{1}V^{*}I_{1}^{*}\frac{yv}{z}-\beta_{2}S^{*}I_{2}^{*}\frac{xw}{u}$
$\displaystyle-\delta_{2}\beta_{2}V^{*}I_{2}^{*}\frac{yw}{u}-(a_{1}+\mu)E_{1}^{*}\frac{z}{v}-(a_{2}+\mu)E_{2}^{*}\frac{u}{w}$
$\displaystyle=$ $\displaystyle\mu S^{*}\Big{(}2-x-\frac{1}{x}\Big{)}+\gamma
V^{*}\Big{(}2-\frac{x}{y}-\frac{y}{x}\Big{)}+\mu
V^{*}\Big{(}3-y-\frac{1}{x}-\frac{x}{y}\Big{)}$
$\displaystyle+\beta_{1}S^{*}I_{1}^{*}\Big{(}3-\frac{1}{x}-\frac{xv}{z}-\frac{z}{v}\Big{)}+\beta_{2}S^{*}I_{2}^{*}\Big{(}3-\frac{1}{x}-\frac{xw}{u}-\frac{u}{w}\Big{)}$
$\displaystyle+\delta_{2}\beta_{2}V^{*}I_{2}^{*}\Big{(}4-\frac{1}{x}-\frac{x}{y}-\frac{yw}{u}-\frac{u}{w}\Big{)}+\delta_{1}\beta_{1}V^{*}I_{1}^{*}\Big{(}4-\frac{1}{x}-\frac{x}{y}-\frac{yv}{z}-\frac{z}{v}\Big{)}.$

Here, clearly $F(x,y,z,u,v,w)\leq 0,$ and the equality holds only for $x=y=1$,
$z=v$, and $u=w$, i.e.,

$\left\\{(x,y,z,u,v,w)\in\Delta:F(x,y,z,u,v,w)=0\right\\}\equiv\left\\{(x,y,z,u,v,w):x=y=1,z=v,u=w\right\\},$

which corresponds to the set

$\Delta^{\prime}=\left\\{(S,V,E_{1},E_{2},I_{1},I_{2}):S=S^{*},V=V^{*},\frac{E_{1}}{E_{1}^{*}}=\frac{I_{1}}{I_{1}^{*}},\frac{E_{2}}{E_{2}^{*}}=\frac{I_{2}}{I_{2}^{*}}\right\\}\subset\Delta.$

It is evident to see that the maximum invariant set on the set
$\Delta^{\prime}$ is the singleton $\left\\{D^{*}\right\\}$, then the endemic
equilibrium $D^{*}$ is globally stable in $\Delta$ by LaSalle’s Invariance
Principle [40]. ∎

###  4.4 Bifurcations

Different dynamical behaviors may occur in a mathematical model for the
variation of the model parameters. The critical parameter value at which
qualitative dynamics change occurs is called a bifurcation point. The
objective of this section is to determine some local bifurcations of the
system (1) with the variation of different parameters.

####  4.4.1 Hopf bifurcation

This section focuses on the local stability and Hopf bifurcation at the
positive equilibrium $D^{*}$ of system (1), which represents the coexistence
of the both strains (native and mutant). To determine the local asymptotic
stability of $D^{*}$, the characteristic equation of the linearized system of
(1) at $D^{*}$ is utilized. The characteristic equation is given by

$C(\lambda)=\lambda^{7}+l_{1}\lambda^{6}+l_{2}\lambda^{5}+l_{3}\lambda^{4}+l_{4}\lambda^{3}+l_{5}\lambda^{2}+l_{6}\lambda+l_{7}=0,$
(17)

where $l_{1},l_{2},l_{3},l_{4},l_{5},l_{6},$ and $l_{7}$ are given in Appendix
A. Now, we define Routh-Hurwitz determinants

$\displaystyle H_{1}=$ $\displaystyle l_{1}\leavevmode\nobreak\ ,\quad
H_{2}=\begin{array}[]{|cc|}l_{1}&l_{3}\\\
1&l_{2}\end{array}\leavevmode\nobreak\ ,\quad
H_{3}=\begin{array}[]{|ccc|}l_{1}&l_{3}&l_{5}\\\ 1&l_{2}&l_{4}\\\
0&l_{1}&l_{3}\end{array}\leavevmode\nobreak\ ,\quad
H_{4}=\begin{array}[]{|cccc|}l_{1}&l_{3}&l_{5}&l_{7}\\\ 1&l_{2}&l_{4}&l_{6}\\\
0&l_{1}&l_{3}&l_{5}\\\ 0&1&l_{2}&l_{4}\end{array}\leavevmode\nobreak\ ,$
$\displaystyle H_{5}=$
$\displaystyle\begin{array}[]{|ccccc|}l_{1}&l_{3}&l_{5}&l_{7}&0\\\
1&l_{2}&l_{4}&l_{6}&0\\\ 0&l_{1}&l_{3}&l_{5}&l_{7}\\\ 0&1&l_{2}&l_{4}&l_{6}\\\
0&0&l_{1}&l_{3}&l_{5}\end{array}\leavevmode\nobreak\ ,\quad
H_{6}=\begin{array}[]{|cccccc|}l_{1}&l_{3}&l_{5}&l_{7}&0&0\\\
1&l_{2}&l_{4}&l_{6}&0&0\\\ 0&l_{1}&l_{3}&l_{5}&l_{7}&0\\\
0&1&l_{2}&l_{4}&l_{6}&0\\\ 0&0&l_{1}&l_{3}&l_{5}&l_{7}\\\
0&0&1&l_{2}&l_{4}&l_{6}\end{array}\leavevmode\nobreak\ ,\quad
H_{7}=l_{7}H_{6}.$

By Routh-Hurwitz criterion, $D^{*}$ is locally asymptotically stable (i.e.
$Re(\lambda)<0$) if and only if
$H_{1}>0,H_{2}>0,H_{3}>0,H_{4}>0,H_{5}>0,H_{6}>0,$ and $H_{7}>0;$ otherwise,
$D^{*}$ becomes unstable.

Further, we determine the occurrence conditions of Hopf bifurcation of system
(1). To study the Hopf bifurcation, bifurcation parameter should be chosen at
first. Among all parameters of system (1), we choose the parameter $\nu_{1}$,
which represents the mutation rate of the native strain. With the other
parameter values given, we can calculate the threshold value of the
bifurcation parameter $\nu_{1}$. In the following, we denote this threshold
value of Hopf bifurcation point as $\nu_{1}=\nu_{1}^{*}$.

By Liu criterion [46], we assume that there is a smooth curve of equilibrium
points $(D(\nu_{1}),\nu_{1}^{*})$ with $D(\nu_{1})=\nu_{1}^{*}$ for system (1)
and $(D^{*},\nu_{1}^{*})$ is a positive equilibrium point. If
$C(\lambda,\nu_{1}^{*})=\lambda^{7}+l_{1}(\nu_{1}^{*})\lambda^{6}+l_{2}(\nu_{1}^{*})\lambda^{5}+l_{3}(\nu_{1}^{*})\lambda^{4}+l_{4}(\nu_{1}^{*})\lambda^{3}+l_{5}(\nu_{1}^{*})\lambda^{2}+l_{6}(\nu_{1}^{*})\lambda+l_{7}(\nu_{1}^{*})$
is the characteristic equation at $(D^{*},\nu_{1}^{*})$, then for a simple
Hopf bifurcation, we have the following conditions:

  * (i)

$l_{7}(\nu_{1}^{*})>0,H_{1}(\nu_{1}^{*})>0,H_{2}(\nu_{1}^{*})>0,H_{3}(\nu_{1}^{*})>0,H_{4}(\nu_{1}^{*})>0,H_{5}(\nu_{1}^{*})>0,$
and $H_{6}(\nu_{1}^{*})=0.$

  * (ii)

$\frac{d}{d\nu_{1}}(H_{6}(\nu_{1}^{*}))\neq 0,$

where
$H_{1}(\nu_{1}^{*}),H_{2}(\nu_{1}^{*}),H_{3}(\nu_{1}^{*}),H_{4}(\nu_{1}^{*}),H_{5}(\nu_{1}^{*}),$
and $H_{6}(\nu_{1}^{*})$ are the Hurwitz determinants at the bifurcation
parameter $\nu_{1}^{*}$.

If the condition (i) holds, then the characteristic polynomial have to meet
the condition for a pair of purely imaginary eigenvalues. Now, to occurrence
the Hopf bifurcation, we need to derive the transversality condition (ii). For
this, we let $\pm i\omega$ be a pair of purely imaginary eigenvalues. Here,
differentiating the characteristic equation (17) with respect to $\nu_{1},$ we
obtain

$\displaystyle(7\lambda^{6}+6l_{1}\lambda^{5}+5l_{2}\lambda^{4}+4l_{3}\lambda^{3}$
$\displaystyle+3l_{4}\lambda^{2}+2l_{5}\lambda+l_{6})\frac{d\lambda}{d\nu_{1}}+\lambda^{6}\frac{dl_{1}}{d\nu_{1}}+\lambda^{5}\frac{dl_{2}}{d\nu_{1}}$
$\displaystyle+\lambda^{4}\frac{dl_{3}}{d\nu_{1}}+\lambda^{3}\frac{dl_{4}}{d\nu_{1}}+\lambda^{2}\frac{dl_{5}}{d\nu_{1}}+\lambda\frac{dl_{6}}{d\nu_{1}}+\frac{dl_{7}}{d\nu_{1}}=0.$

Further, we obtain

$\displaystyle\Bigg{(}\frac{d\lambda}{d\nu_{1}}\Bigg{)}^{-1}=-\frac{7\lambda^{6}+6l_{1}\lambda^{5}+5l_{2}\lambda^{4}+4l_{3}\lambda^{3}+3l_{4}\lambda^{2}+2l_{5}\lambda+l_{6}}{\lambda^{6}\frac{dl_{1}}{d\nu_{1}}+\lambda^{5}\frac{dl_{2}}{d\nu_{1}}+\lambda^{4}\frac{dl_{3}}{d\nu_{1}}+\lambda^{3}\frac{dl_{4}}{d\nu_{1}}+\lambda^{2}\frac{dl_{5}}{d\nu_{1}}+\lambda\frac{dl_{6}}{d\nu_{1}}+\frac{dl_{7}}{d\nu_{1}}}.$

Furthermore, we have

$\displaystyle
sign\Bigg{[}\frac{d(Re(\lambda))}{d\nu_{1}}\Bigg{]}_{\lambda=i\omega,H_{6}=0}$
$\displaystyle=sign\Bigg{[}Re\Bigg{(}\frac{d\lambda}{d\nu_{1}}\Bigg{)}^{-1}\Bigg{]}_{\lambda=i\omega,H_{6}=0}$
$\displaystyle=sign[\Upsilon],$

where,

$\displaystyle\Upsilon$
$\displaystyle=Re\Bigg{[}\frac{(7\omega^{6}-5l_{2}\omega^{4}+3l_{4}\omega^{2}-l_{6})+i(-6l_{4}\omega^{5}+4l_{3}\omega^{2}-2l_{5}\omega)}{(-\omega^{6}\frac{dl_{1}}{d\nu_{1}}+\omega^{4}\frac{dl_{3}}{d\nu_{1}}-\omega^{2}\frac{dl_{5}}{d\nu_{1}}+\frac{dl_{7}}{d\nu_{1}})+i(\omega^{5}\frac{dl_{2}}{d\nu_{1}}-\omega^{3}\frac{dl_{4}}{d\nu_{1}}+\omega\frac{dl_{6}}{d\nu_{1}})}\Bigg{]}=\frac{M_{1}M_{3}+M_{2}M_{4}}{M_{3}^{2}+M_{4}^{2}},$
$\displaystyle M_{1}$
$\displaystyle=7\omega^{6}-5l_{2}\omega^{4}+3l_{4}\omega^{2}-l_{6},\quad
M_{2}=-6l_{4}\omega^{5}+4l_{3}\omega^{2}-2l_{5}\omega,$ $\displaystyle M_{3}$
$\displaystyle=-\omega^{6}\frac{dl_{1}}{d\nu_{1}}+\omega^{4}\frac{dl_{3}}{d\nu_{1}}-\omega^{2}\frac{dl_{5}}{d\nu_{1}}+\frac{dl_{7}}{d\nu_{1}},\quad
M_{4}=\omega^{5}\frac{dl_{2}}{d\nu_{1}}-\omega^{3}\frac{dl_{4}}{d\nu_{1}}+\omega\frac{dl_{6}}{d\nu_{1}}.$

If $M_{1}M_{3}+M_{2}M_{4}>0,$ then
$sign\Big{[}\frac{d(Re(\lambda))}{d\nu_{1}}\Big{]}_{\nu_{1}=\nu_{1}^{*}}>0$
and the transversality condition (ii) holds. Summarizing the above discussion,
we obtain the subsequent theorem:

######  Theorem 4.6.

For the existing positive equilibrium $D^{*}$ of system (1), the system (1)
around $D^{*}$ enters into Hopf bifurcation when $\nu_{1}$ crosses through
$\nu_{1}^{*}$.

####  4.4.2 Transcritical bifurcation

We see that Eq. (6) has a zero eigenvalue when either $R_{1}=1$ or $R_{2}=1$.
Thus, the system (1) may undergo a transcritical bifurcation at $D^{0}$ when
either $R_{1}=1$ or $R_{2}=1$. In this subsection, we establish conditions on
the parameters using Theorem 4.1 from Castillo-Chavez and Song [16] and center
manifold theory [26]. For the transcritical bifurcation, we establish the
following theorem:

######  Theorem 4.7.

  * 1.

Assume $R_{1}<1$, the system (1) undergoes a transcritical bifurcation near
$D^{0}$, when $R_{2}=1$.

  * 2.

Assume $R_{2}<1$, the system (1) undergoes a transcritical bifurcation near
$D^{0}$, when $R_{1}=1$.

###### Proof.

1\. We choose $\beta_{2}$ as a bifurcation parameter. By solving $R_{2}=1,$ we
obtain

$\beta_{2}=\beta_{2}^{*}=\frac{\mu\left(a_{2}+\mu\right)\left(\alpha_{2}+d_{2}+\mu\right)(\gamma+\mu+p)}{a_{2}\Lambda\left(\gamma+\mu+\delta_{2}p\right)}.$

It can easily be obtained that the Jacobian $J_{(D^{0},\beta_{2}^{*})}$
evaluated at $D^{0}$ and $\beta_{2}=\beta_{2}^{*}$ has a simple zero
eigenvalue and other eigenvalues have negative sign. Hence $D^{0}$ is a non-
hyperbolic equilibrium, when $\beta_{2}=\beta_{2}^{*}.$ Now, we calculate a
right eigenvector $W=(w_{1},w_{2},w_{3},w_{4},w_{5},w_{6})$ and a left
eigenvector $V=(v_{1},v_{2},v_{3},v_{4},v_{5},v_{6})$ associated to the zero
eigenvalues. Here

$\displaystyle w_{1}$
$\displaystyle=-\frac{\left(a_{2}+\mu\right)\left(\alpha_{2}+d_{2}+\mu\right)\left((\gamma+\mu)^{2}+\gamma\delta_{2}p\right)}{a_{2}\mu(\gamma+\mu+p)\left(\gamma+\mu+\delta_{2}p\right)},$
$\displaystyle w_{2}$
$\displaystyle=-\frac{p\left(a_{2}+\mu\right)\left(\alpha_{2}+d_{2}+\mu\right)\left(\gamma+\mu+\delta_{2}(\mu+p)\right)}{a_{2}\mu(\gamma+\mu+p)\left(\gamma+\mu+\delta_{2}p\right)},$
$\displaystyle w_{3}$ $\displaystyle=0,\quad
w_{4}=\frac{\mu+d_{2}+\alpha_{2}}{a_{2}},\quad w_{5}=0,\quad w_{6}=1,\quad
w_{7}=0,$ $\displaystyle v_{1}$ $\displaystyle=0,\quad v_{2}=0,\quad
v_{3}=\frac{a_{1}\nu_{1}}{(\mu+a_{1})(\mu+d_{1}+\alpha_{1}+\nu_{1})(1-R_{1})},$
$\displaystyle v_{4}$ $\displaystyle=1,\quad
v_{5}=\frac{\nu_{1}}{(\mu+d_{1}+\alpha_{1}+\nu_{1})(1-R_{1})},\quad
v_{6}=1,\quad v_{7}=0.$

Now from Theorem 4.1 of [16], we need to calculate the bifurcation constants
$a$ and $b$. For system (1), the associated non-zero partial derivatives of
$f$ (evaluated at $D^{0}$,
$x_{1}=S,x_{2}=V,x_{3}=E_{1},x_{4}=E_{2},x_{5}=I_{1},x_{6}=I_{2},x_{7}=R_{1}$)
are given by

$\displaystyle a$
$\displaystyle=2v_{3}w_{1}w_{5}\frac{\partial^{2}f_{3}}{\partial S\partial
I_{1}}+2v_{3}w_{2}w_{5}\frac{\partial^{2}f_{3}}{\partial V\partial
I_{1}}+2v_{4}w_{1}w_{6}\frac{\partial^{2}f_{4}}{\partial S\partial
I_{2}}+2v_{4}w_{2}w_{6}\frac{\partial^{2}f_{4}}{\partial V\partial
I_{2}}+2v_{4}w_{6}w_{7}\frac{\partial^{2}f_{4}}{\partial I_{2}\partial R_{1}}$
$\displaystyle=2v_{4}w_{6}(w_{1}+\delta_{2}w_{2})\beta_{2}^{*}<0,$
$\displaystyle b$
$\displaystyle=2v_{4}w_{6}\frac{\partial^{2}f_{4}}{\partial\beta_{2}\partial
I_{2}}=2v_{4}w_{6}(S^{0}+\delta_{2}V^{0})>0.$

Since the coefficient $a$ is negative and $b$ is positive, the direction of
the bifurcation of system (1) at $\beta_{2}=\beta_{2}^{*}$ is forward.
2\. By choosing $\beta_{1}$ as a bifurcation parameter. By solving $R_{1}=1,$
we obtain

$\beta_{1}=\beta_{1}^{*}=\frac{\mu\left(a_{1}+\mu\right)(\gamma+\mu+p)\left(\alpha_{1}+d_{1}+\mu+\nu_{1}\right)}{a_{1}\Lambda\left(\gamma+\mu+\delta_{1}p\right)}.$

Following simple procedure of previous case, we obtain that $D^{0}$ is a non-
hyperbolic equilibrium, when $\beta_{1}=\beta_{1}^{*}.$ Now, we calculate a
right eigenvector $W_{1}=(w_{11},w_{22},w_{33},w_{44},w_{55},w_{66})$ and a
left eigenvector $V_{1}=(v_{11},v_{22},v_{33},v_{44},v_{55},v_{66})$
associated to the zero eigenvalues. Here

$\displaystyle w_{11}$
$\displaystyle=-\frac{\Lambda\beta_{2}((\gamma+\mu)^{2}+p\gamma\delta_{2})\nu_{1}}{\alpha_{1}\mu(p+\gamma+\mu)^{2}(\mu+d_{2}+\alpha_{2})(\mu+a_{2})(1-R_{2})}-\frac{\left(a_{1}+\mu\right)\left(\alpha_{1}+d_{1}+\mu+\nu_{1}\right)\left((\gamma+\mu)^{2}+\gamma\delta_{1}p\right)}{a_{1}\alpha_{1}(\gamma+\mu+p)\left(\gamma+\mu+\delta_{1}p\right)},$
$\displaystyle w_{22}$
$\displaystyle=-\frac{p(\gamma+\mu+(p+\mu)\delta_{2})}{(p+\gamma+\mu)\alpha_{1}}\Big{(}\frac{\Lambda\beta_{2}\nu_{1}}{\mu(p+\gamma+\mu)(\mu+d_{2}+\alpha_{2})(1-R_{2})}+\frac{(\mu+a_{1})(\mu+d_{1}+\alpha_{1}+\nu_{1})}{a_{1}(\gamma+\mu+p\delta_{2})}\Big{)},$
$\displaystyle w_{33}$
$\displaystyle=\frac{\mu\left(\alpha_{1}+d_{1}+\mu+\nu_{1}\right)}{a_{1}\alpha_{1}},\quad
w_{44}=\frac{\Lambda\beta_{2}(\gamma+\mu+p\delta_{2})\nu_{1}}{\alpha_{1}(p+\gamma+\mu)(\mu+d_{2}+\alpha_{2})(\mu+a_{2})(1-R_{2})},\quad
w_{55}=\frac{\mu}{\alpha_{1}},$ $\displaystyle w_{66}$
$\displaystyle=\frac{\mu\nu_{1}}{\alpha_{1}(\mu+d_{2}+\alpha_{2})(1-R_{2})},\quad
w_{77}=1\quad v_{11}=0,\quad v_{22}=0,\quad
v_{33}=\frac{a_{1}}{a_{1}+\mu},\quad v_{44}=0,$ $\displaystyle v_{55}$
$\displaystyle=1,\quad v_{66}=0,\quad v_{77}=0.$

Similarly, as in previous case, we have

$\displaystyle a$
$\displaystyle=2v_{33}w_{11}w_{55}\frac{\partial^{2}f_{3}}{\partial S\partial
I_{1}}+2v_{33}w_{22}w_{55}\frac{\partial^{2}f_{3}}{\partial V\partial I_{1}}$
$\displaystyle=2v_{33}w_{11}w_{55}\beta_{1}^{*}+2v_{33}w_{22}w_{55}\delta_{1}\beta_{1}^{*},$
$\displaystyle=2v_{33}w_{55}\beta_{1}^{*}(w_{11}+\delta_{1}w_{22})<0,$
$\displaystyle b$
$\displaystyle=2v_{33}w_{55}\frac{\partial^{2}f_{3}}{\partial\beta_{1}\partial
I_{1}}=2v_{33}w_{55}(S^{0}+\delta_{1}V^{0})>0.$

Therefore again in this case, the direction of the bifurcation of system (1)
at $\beta_{1}=\beta_{1}^{*}$ is forward. ∎

##  5 Implications for disease control

For vaccine-preventable diseases, not all susceptible individuals could be
immunized for numerous reasons, such as they are too young to be vaccinated
(vaccinating young children or infants can, sometimes, harm the
children/infants), they have weakened immune system and co-morbidities (and
vaccinating may make their prognosis worse), they are of higher age, or they
opt out for religious or traditional reasons, cultural. However, the question
is what is the smallest proportion of those individuals we can vaccinate that
we need to vaccinate so that those we cannot vaccinate can be protected from
emerging severe disease or dying of the disease. The idea of herd immunity in
the disease dynamics is related to the indirect protection against acquiring
of infectious disease, which members of the community obtain when a large
percentage of the populace has become protected to the contagious disease due
to natural recovery from prior infection or vaccination [3, 4, 21]. The
outcome of herd immunity is that persons who are not immune (e.g., those who
have not been infected yet or cannot be vaccinated) obtain some defense
against acquiring the infection. The fastest and safest and way to attain herd
immunity is vaccination. It should, however, be stated that Sweden implemented
the other procedure for achieving herd immunity in the COVID-19 dynamics in
Sweden [23]. In other words, the Swedish public health agencies intended to
achieve herd immunity by not employing the basic society transmission
reduction strategies (e.g., community lockdowns, social distancing, contact
tracing use of face masks in public, etc.) employed in almost every nation or
community that is hard-hit with the COVID-19 pandemic, opting, instead, to let
persons to attain disease and, hopefully, recover from it. In this section, a
theoretical condition for achieving community-wide vaccine-induced herd
immunity is obtained. Theorem 4.2 has significant public health implications.
It reveals that if the imperfect vaccine has sufficient efficacy and coverage
rate to make $R_{0}<1,$ COVID-19 will be eradicated from society. The global
stability of the disease-free equilibrium (Theorem 4.3) for $R_{0}<1$ confirms
that such epidemics do not hit. This means $R_{0}$ is an appropriate
combination of parameters to measure the efficiency of a vaccination campaign.

###  5.1 Herd immunity

Not every person in a given population expects to be immunized in order to
eradicate the disease. A fraction of people with immunity in the given
population required to stop an epidemic is named herd immunity. Let $\rho$
represent the fraction of the vaccinated population at $D^{0}$ (the disease-
free equilibrium). Then,

$\rho=\frac{p}{\mu+\gamma+p}.$

In the absence of vaccination, i.e., when $p=0,$ the basic reproduction number
is given by $R_{0wv}$. Hence, we can write

$R_{0}=\max\left\\{R_{1wv}(1-(1-\delta_{1})\rho),R_{2wv}(1-(1-\delta_{2})\rho)\right\\},$
(18)

with $R_{1}=R_{1wv}(1-(1-\delta_{1})\rho),$ and
$R_{2}=R_{2wv}(1-(1-\delta_{2})\rho).$ It is noted that $R_{1}\leq
R_{1wv},\leavevmode\nobreak\ R_{2}\leq R_{2wv},$ and thus $R_{0}\leq R_{0wv}.$
The equality holds only $\rho=0$ (i.e., $p=0$) or $\delta_{1}=\delta_{2}=1.$
This indicates that the vaccine, even not $100\%$ effective, will certainly
reduce the basic reproduction number of the disease. Since $R_{0}<1$ is a
necessary and sufficient condition for the eradication of disease (Theorems
4.2 and 4.3), hence it follows from (18) that

$\rho>\max\left\\{\frac{1}{1-\delta_{1}}\Bigg{(}1-\frac{1}{R_{1wv}}\Bigg{)},\frac{1}{1-\delta_{2}}\Bigg{(}1-\frac{1}{R_{2wv}}\Bigg{)}\right\\}=\rho_{critical}$
(19)

is also a necessary and sufficient condition for disease elimination. Here,
$\rho_{critical}$ signifies herd immunity. Although, this outcome could be
achieved in the case of continuous vaccination and that herd immunity is
attained if the vaccination rate is large enough such that $\rho$, the
fraction of vaccinated individuals at the disease-free equilibrium, exceeds
the critical value $\rho_{critical}$. Combining Theorems 4.2 and 4.3, we
obtain the following consequence:

######  Proposition 5.1.

COVID-19 can be eradicated from the population if $\rho>\rho_{critical}.$

The inequality (19) can be expressed in terms of the vaccination rate $p$.
This is done by noting, first of all, that $R_{0}$ is a decreasing function of
$p$,

$\frac{dR_{0}}{dp}=\max\left\\{\frac{-\Lambda(\gamma+\mu)(1-\delta_{1})a_{1}\beta_{1}}{\mu(p+\gamma+\mu)^{2}(\mu+a_{1})(\mu+d_{1}+\alpha_{1}+\nu_{1})},\frac{-\Lambda(\gamma+\mu)(1-\delta_{2})a_{2}\beta_{2}}{\mu(p+\gamma+\mu)^{2}(\mu+a_{2})(\mu+d_{2}+\alpha_{2})}\right\\}<0$

and so it is minimized by letting $p$ go to infinity. Taking the limit as $p$
approaches infinity, we see that this expression is always greater than
$\max\left\\{\delta_{1}R_{1wv},\delta_{2}R_{2wv}\right\\}$. Thus, if
$\max\left\\{\delta_{1}R_{1wv},\delta_{2}R_{2wv}\right\\}>1,$ then no amount
of vaccination can make $R_{0}$ smaller than unity. Alternatively, if
$\max\left\\{\delta_{1}R_{1wv},\delta_{2}R_{2wv}\right\\}<1$, then the
condition

$p>\max\left\\{\frac{(\gamma+\mu)(R_{1wv}-1)}{1-\delta_{1}R_{1wv}},\frac{(\gamma+\mu)(R_{2wv}-1)}{1-\delta_{2}R_{2wv}}\right\\}=p_{critical}$
(20)

gives $R_{0}<1$. Of course, this condition assumes $R_{0}>1$, since disease
elimination follows without vaccination if $R_{0}<1$ (by Theorems 4.2 and 4.3
and the fact that $R_{0}\leq R_{0wv}$). It is easy to show that from (20), we
obtain $R_{0}<1$ if $p>p_{critical}$, and $R_{0}>1$ if $p<p_{critical}$. Thus,
we have established the following result:

######  Proposition 5.2.

If $\max\left\\{\delta_{1}R_{1wv},\delta_{2}R_{2wv}\right\\}<1$ and
$p>p_{critical}$, then COVID-19 will be eliminated from the community. If
$\max\left\\{\delta_{1}R_{1wv},\delta_{2}R_{2wv}\right\\}>1$, then no amount
of vaccination would be able to prevent the COVID-19 outbreak in the
community.

##  6 Numerical illustration

In this section, we investigate the dynamics of system (1) numerically for
different sets of parameter values. Such investigations aim to determine the
effect of varying the values of the different parameters and support the
obtained theoretical results. It is observed that the hypothetical values of
parameters given in Table 2 are biologically feasible. However, to verify the
bifurcations and different dynamical behavior of system (1), some parameters
are varied differently from Table 2.

Table 2: Numerical values of parameters. Parameters |  Value with unit |  Reference
---|---|---
$\Lambda$ |  10000 people per week |  Assumed
$\beta_{1}$ |  0.000001 |  Assumed
$\beta_{2}$ |  0.0000003 |  Assumed
$p$ |  0.02 week-1 |  Assumed
$\mu$ |  0.0003 week-1 |  $1/(65\times 48)$
$1-\delta_{1}$ |  0.75 (dimensionless) |  Assumed that vaccine efficacy is 75% to native strain
$1-\delta_{2}$ |  0.40 (dimensionless) |  Assumed that vaccine efficacy is 40% to mutant strain
$\gamma$ |  1/32 week-1 |  Assumed that loss of immunity of vaccinated individuals after 32 weeks
$a_{1}$ |  1 week-1 |  Incubation (1 Week) [41]
$a_{2}$ |  1 week-1 |  Incubation (1 Week) [41]
$\alpha_{1}$ |  1/2 week-1 |  Recovery (2 Weeks) [34]
$\alpha_{2}$ |  1/2 week-1 |  Recovery (2 Weeks) [34]
$d_{1}$ |  0.0006 week-1 |  Assumed
$d_{2}$ |  0.0006 week-1 |  Assumed
$1-\delta_{3}$ |  0.90 (dimensionless) |  Assumed that recovered individuals have 90% immunity against mutant strain
$\nu_{1}$ |  0.3 week-1 |  Assumed

The dynamics of system (1) are simulated using MATLAB 2018a. Figure 2(a) shows
that for $R_{0}=\max\left\\{R_{1},R_{2}\right\\}=0.9735<1,$ DFE $(D^{0})$.is
stable. Figure 2(b) represents that the mutant dominant equilibrium ($D^{2}$)
is stable for $R_{1}<1$ and $R_{2}>1$. Further, Figure 3 shows that mutation
rate ($\nu_{1}$) can change the dynamics of positive equilibrium of system
(1). Figure 3(a) represents that the positive equilibrium is stable. Figure
3(b) and (c) ensure that system (1) losses its stability around the positive
equilibrium for high mutation rate and periodic solution occurs. Furthermore,
Figure 3(d) shows that the positive equilibrium regains its stability for
further higher mutation rate. This type of dynamics is called endemic bubble
[45].

(a) (b)

Figure 2: (a) The solution $I_{1}$ and $I_{2}$ of system (1) for the
parametric values $\beta_{1}=0.000000033,\beta_{2}=0.000000017,$ and others
from Table 2. (b) The solution $I_{1}$ and $I_{2}$ of system (1) for
$\beta_{1}=0.0000000033,\beta_{2}=0.0000003,$ and other parametric values from
Table 2.

(a) (b) (c) (d)

Figure 3: (a) The solution $I_{1}$ and $I_{2}$ of system (1) for the
parametric values in Table 2 and $\nu_{1}=0.3$. (b) The solution $I_{1}$ and
$I_{2}$ of system (1) for $\nu_{1}=0.5$. (c) The solution $I_{1}$ and $I_{2}$
of system (1) for $\nu_{1}=0.7$. (d) The solution $I_{1}$ and $I_{2}$ of
system (1) for $\nu_{1}=0.90$.

Furthermore, Figure 4 shows the positive equilibrium is locally asymptotically
stable for small and larger enough values of $\nu_{1}$, but unstable for
intermediate values of $\nu_{1}$. In this way, we get the bifurcation diagram
in Figure 4, which we call an endemic bubble. As shown in the bifurcation
diagram (Figure 4), for a lower range of values of $\nu_{1}$, both strains
persist in the environment and system (1) is asymptotically stable. For a
range of $\nu_{1}$, the periodic oscillations (limit cycle) will appear, but
for higher values of $\nu_{1},$ system (1) regains its stability. Figure 5
illustrates the Hopf bifurcation diagram with respect to the parameter
$\delta_{3}$. It shows that the system (1) is locally asymptotically stable
for lower values of $\delta_{3}$ and periodic solutions (Hopf bifurcation)
appear for higher values of $\delta_{3}$.

(a) (b)

Figure 4: Bifurcation diagram (endemic bubble) with respect to mutation rate
($\nu_{1}$), other parametric values remain same as in Table 2. The blue color
shows the upper limit of the limit cycle and red color shows the lower limit
of the limit cycle.

(a) (b)

Figure 5: Hopf bifurcation diagram with respect to $\delta_{3}$. We keep
$\nu_{1}=0.4$ and other parametric values same as in Table 2. The blue color
shows the upper limit of the limit cycle and red color shows the lower limit
of the limit cycle.

##  7 Case study on COVID-19 data in India

Here we estimate the unknown parameters of the system (1) on the cumulative
cases, and cumulative deaths of COVID-19 in India from March 1, 2021 to
September 27, 2021 by using the MCMC algorithm [29, 30, 1]. We collect the
data of cumulative cases and cumulative deaths of COVID-19 for March 1, 2021
to September 27, 2021 (31 weeks) from the WHO website [70]. By estimating the
parameters, we estimate the mean values, standard deviation, and Geweke values
of some parameters of the system (1). The cumulative cases can be given as

$\frac{dC}{dt}=a_{1}E_{1}+a_{2}E_{2},$ (21)

where $C(t)$ represents the cumulative cases, and the cumulative deaths can be
given as

$\frac{dD}{dt}=d_{1}I_{1}+d_{2}I_{2},$ (22)

where $D(t)$ represents the cumulative deaths.

We use MCMC method for 20000 simulations to fit the Eq. (21) and estimate the
parameters. Figure 6 represents a good fitting between the cumulative reported
cases of COVID-19 and the model solution, well suggesting the epidemic trend
in India. Some of the parameters of the system (1) are taken either from the
literature or presumed based on publicly-available COVID-19 associated
information. We assume parametric values of
$a_{1},a_{2},\mu,\delta_{1},\delta_{2},$ and $\delta_{3}$ same as given in
Table 2 and initial conditions given in Table 4.

(a) (b)

Figure 6: Fitting result of the cumulative cases and cumulative deaths of
COVID-19. (a) The blue solid boxes represent the actual reported cumulative
cases and the orange curve represents the model output. (b) The blue solid
boxes represent the actual reported cumulative deaths and the orange curve
represents the model output. For the different colors in the Figure, refer to
the web version of the paper.

By using MCMC method, we acquire the values of the parameters
$\beta_{1},\beta_{2},\alpha_{1},\alpha_{2},\nu_{1},d_{1},$ and $d_{2}$ with
MCMC chain of the time evolution of the cumulative cases and comparison with
the confirmed cases of COVID-19 in India. We compute the mean values, standard
deviation, and Geweke values of these parameters (refer in Table 3).

Table 3: Estimated values of parameters by MCMC method. Parameters |  Mean value |  Standard deviation |  Geweke value |  Reference
---|---|---|---|---
$\Lambda$ |  320000 |  - |  - |  Assumed
$\beta_{1}$ |  2.999e-07 |  5.3206e-13 |  0.99 |  MCMC
$\beta_{2}$ |  3.7618e-08 |  3.145e-14 |  0.95 |  MCMC
$p$ |  0.00002 |  - |  - |  Assumed
$\gamma$ |  0.01 |  - |  - |  Assumed
$\alpha_{1}$ |  0.45 |  5.5791e-03 |  0.99 |  MCMC
$\alpha_{2}$ |  0.2 |  1.6627e-03 |  0.98 |  MCMC
$\nu_{1}$ |  0.03 |  3.265e-04 |  0.99 |  MCMC
$d_{1}$ |  0.007 |  5.347e-05 |  0.99 |  MCMC
$d_{2}$ |  0.00003 |  3.016e-7 |  0.99 |  MCMC
Table 4: Initial conditions for system (1) with respect to COVID-19 in India. S(0) |  V(0) |  $E_{1}(0)$ |  $E_{2}(0)$ |  $I_{1}(0)$ |  $I_{2}(0)$ |  $R_{1}(0)$ |  $R_{2}(0)$
---|---|---|---|---|---|---|---
$23\times 10^{5}$ |  12256337 |  $10\times 10^{4}$ |  $60\times 10^{4}$ |  $10\times 10^{3}$ |  $80\times 10^{3}$ |  $60\times 10^{2}$ |  $20\times 10^{3}$

###  7.1 The effects of different efficacies of vaccine and mutation rate on
cumulative cases and cumulative deaths over time

The system (1) is simulated to assess the population-level impact of the
imperfect anti-COVID-19 vaccine in India. The population-level impact of the
vaccine efficacy to both strains (native and mutant) on the burden of the
pandemic is examined first of all. Firstly, we consider different values of
the vaccine efficacy to the native strain. The system (1) is then simulated
using the baseline parameter values in Table 2, 3, and different values of the
vaccine efficacy $(1-\delta_{1})$ to the native strain. The results obtained,
depicted in Figure 7(a) show that, for the vaccine efficacy
$(1-\delta_{1})=75\%$ (assumed), 33,652,745 cumulative cases (the red curve)
has reported by September 27, 2021. Predictions show that the cumulative cases
would be recorded 40,350,000, by February 7, 2022 (31 weeks after September
27, 2021). The simulations, further, show a reduction with increasing values
of the vaccine efficacy $(1-\delta_{1})$ from its baseline value. In
particular, if we consider the vaccine efficacy as $(1-\delta_{1})=90\%,$ then
39,030,000 cumulative cases would be recorded by February 7, 2022,
representing only 3.27% reduction. Figure 7(b) represents that, for the
vaccine efficacy $(1-\delta_{2})=40\%$ (assumed) to the mutant strain,
33,652,745 cumulative cases (the yellow curve) has reported by September 27,
2021. Simulations show that the cumulative cases would be recorded 41,520,000,
by February 7, 2022. These simulations show that a reduction in the cumulative
cases with increasing values of the vaccine efficacy $(1-\delta_{2})$.
Particularly, if we increase the vaccine efficacy $(1-\delta_{2})=75\%$, then
41,220,000 cumulative cases would be recorded by February 7, 2022. This
represents 0.7% reduction in the cumulative cases. Furthermore, by increasing
the vaccine efficacy $(1-\delta_{2})=90\%$, the cumulative cases would be
decreased 0.9%.

The system (1) is simulated for different values of $1-\delta_{3}$ (immunity
against mutant strain). Figure 7(c) represents that, if recovered individuals
by native strain have 90% immunity against mutant strain
($(1-\delta_{3})=90\%$) (assumed) to the mutant strain, 33,652,745 cumulative
cases (the blue curve) has reported by September 27, 2021. Simulations show
that the cumulative cases would be recorded 40,350,000, by February 7, 2022.
These simulations show that an increase in the cumulative cases with
decreasing value of $1-\delta_{3}$. Particularly, if we consider
$(1-\delta_{3})=75\%$, then 40,900,000 cumulative cases would be recorded by
February 7, 2022. This represents 1.36% increase in the cumulative cases. If
we consider $(1-\delta_{3})=50\%$, then 42,920,000 cumulative cases would be
recorded by February 7, 2022, representing 6.36% increase in the cumulative
cases. Furthermore, if $(1-\delta_{3})=25\%$, then 48,140,000 cumulative cases
would be recorded by February 7, 2022, representing 19.3% increase in the
cumulative cases.

Furthermore, simulations were also carried out to assess the impact of
mutation rate on the disease dynamics of COVID-19. Figure 7(d) represents
that, for the baseline value of the mutation rate $(\nu_{1}=0.03)$, 33,652,745
cumulative cases (the red curve) has reported by September 27, 2021.
Simulations show that the cumulative cases would be 40,370,000 by February 7,
2022. These simulations show an increase in the cumulative cases with
increasing values of $\nu_{1}$. Particularly, an increase in $\nu_{1}=0.05$,
then 40,380,000 cumulative cases would be recorded by February 7, 2022. This
represents a 0.02% increase in the cumulative cases. Furthermore, an increase
in $\nu_{1}=0.09$, then 40,480,000 cumulative cases would be recorded by
February 7, 2022, representing 0.27% increase in the cumulative cases.

(a) (b) (c) (d)

Figure 7: Assessment of the impacts of different parameters on the cumulative
cases of COVID-19 in India. Simulations of the system (1) shows the cumulative
cases of COVID-19 in India, as a function of time, (a) for different values of
vaccine efficacy to native strain ($1-\delta_{1}$). (b) for different values
of vaccine efficacy to mutant ($1-\delta_{2}$). (c) for different values of
immunity against mutant strain ($1-\delta_{3}$). (d) for different values of
mutation rate ($\nu_{1}$). The baseline parameter values are used from Tables
2 and 3.

To examine the dependence of the end time (February 7, 2022) of the epidemic
on the vaccine efficacies to both native and mutant strains, i.e.
$(1-\delta_{1})$, $(1-\delta_{2})$, mutation rate $(\nu_{1})$, and immunity
against mutant strain $(1-\delta_{3})$, we sketch the contour plots of the
total number of the cumulative cases of COVID-19, with respect to
$(1-\delta_{1})$ and $(1-\delta_{2})$ in Figure 8(a); $(1-\delta_{3})$ and
$\nu_{1}$ in Figure 8(b); $(1-\delta_{2})$ and $\nu_{1}$ in Figure 8(c);
$(1-\delta_{2})$ and $(1-\delta_{3})$ in Figure 8(d); respectively. In order
to assess the combined effect of the parameters, the other parameters remain
fixed when we vary two parameters. The results indicate that increasing the
vaccine efficacies $(1-\delta_{1})$ and $(1-\delta_{2})$, the cumulative cases
would be reduced (Figure 8(a)). It can also be observed that the vaccine
efficacy to the native strain $(1-\delta_{1})$ is more influential than the
vaccine efficacy to the mutant strain $(1-\delta_{2})$, in controlling the
total number of cases. In the same way, the result in Figure 8(b) shows that
the cumulative number of cases would be reduced for the higher immunity
against mutant strain $(1-\delta_{3})$ and lower value of mutation rate
$(\nu_{1})$. The result also shows that cumulative cases would be higher for
the lower value of $(1-\delta_{3})$ and lower value of $(\nu_{1})$. The result
in Figure 8(c) shows that the cumulative number of cases would be reduced for
the higher vaccine efficacy to the mutant strain $(1-\delta_{2})$ and lower
value of mutation rate $(\nu_{1})$. The result also shows that cumulative
cases would be higher for the lower value of $(1-\delta_{2})$ and higher value
of $(\nu_{1})$. Figure 8(d) represents that the cumulative number of cases
would be reduced only for higher immunity against mutant strain
$(1-\delta_{3})$. This result makes it wonder that cumulative cases would be
higher even for the higher value of the vaccine efficacy to the mutant strain
$(1-\delta_{2})$ if there is a low value of immunity against mutant strain.
This means there may exist an appropriate combination of these parameters to
ensure fewer confirmed cases. Thus, it would be interesting to consider an
optimal strategy for supplying vaccines to minimize the cumulative number of
cases.

(a) (b) (c) (d)

Figure 8: Contour plots of the cumulative cases, as a function of different
parameters: (a) as a function of vaccine efficacy ($1-\delta_{1}$) to native
strain and vaccine efficacy ($1-\delta_{2}$) to mutant strain; (b) as a
function of immunity against mutant strain ($1-\delta_{3}$) and mutation rate
($\nu_{1}$); (c) as a function of vaccine efficacy ($1-\delta_{2}$) to mutant
strain and mutation rate ($\nu_{1}$); (d) as a function of vaccine efficacy
($1-\delta_{2}$) to mutant strain and immunity against mutant strain
($1-\delta_{3}$). The baseline parameter values are used from Tables 2 and 3.

Furthermore, the impact of vaccine efficacies on cumulative deaths is
examined. The result in Figure 9(a) shows that for the vaccine efficacy
$(1-\delta_{1})=75\%$ (assumed), 446918 cumulative deaths (red curve) has
reported by September 27, 2021. Predictions show that the cumulative deaths
would be recorded at 534431 by February 7, 2022 (31 weeks after September 27,
2021). The simulations further show a reduction with increasing values of the
vaccine efficacy $(1-\delta_{1})$ from its baseline value. In particular, if
we consider the vaccine efficacy as $(1-\delta_{1})=90\%,$ then 506562
cumulative deaths would be recorded by February 7, 2022, representing only a
5.2% reduction. Figure 9(b) represents that, for the vaccine efficacy
$(1-\delta_{2})=40\%$ (assumed) to the mutant strain, 446918 cumulative deaths
(yellow curve) has reported by September 27, 2021. Simulations show that the
cumulative deaths would be recorded at 543330 by February 7, 2022. Further
simulations show a reduction in the cumulative deaths with increasing values
of the vaccine efficacy $(1-\delta_{2})$. Particularly, if we increase the
vaccine efficacy $(1-\delta_{2})=75\%$, then 537342 cumulative deaths would be
recorded by February 7, 2022, representing a 1.1% reduction in the cumulative
deaths. Furthermore, by increasing the vaccine efficacy $(1-\delta_{2})=90\%$,
529420 cumulative deaths would be recorded, representing 2.5% reduction.

(a) (b)

Figure 9: Assessment of the impacts of different parameters on the cumulative
deaths of COVID-19 in India. Simulations of the system (1) show the cumulative
deaths of COVID-19 in India, as a function of time, (a) for different values
of vaccine efficacy to native strain ($1-\delta_{1}$); (b) for different
values of vaccine efficacy to mutant ($1-\delta_{2}$). The baseline parameter
values are used from Tables 2 and 3.

###  7.2 Impact of mutation rate on the dynamics of strains and infected
population over time

The impact of mutation rate ($\nu_{1}$) has been analyzed on the dynamics of
strains and total infected population for COVID-19. We can explicitly explore
how the numbers of infected individuals $I_{1}$ and $I_{2}$ depend on the
mutation parameter $\nu_{1}$. As shown in the bifurcation diagram (Figure 10),
for the range of mutation rate, $0<\nu_{1}\leq 1.7$, both strains persist in
the environment, and disease is asymptotically stable. For $1.7<\nu_{1}\leq
3.3$, the periodic oscillations (limit cycle) will appear, i.e. the disease
outbreak will occur repeatedly. However, for a higher value of $\nu_{1},$ i.e.
for $\nu_{1}>3.3,$ the periodic solutions disappear and the disease again
becomes stable. This dynamical phenomenon has been illustrated in Figure 11,
for the values $\nu_{1}=1,2.5$ and $3.53$. Furthermore, we have plotted the
total infected population $(E_{1}+E_{2}+I_{1}+I_{2})$ over time for different
values of $\nu_{1}$ in Figure 12. This result gives a wondering dynamics over
a long time. We can easily observe that if the mutation rate $(\nu_{1})$
increases from 0.01 to 0.09, then the infected population persists at a lower
level (blue and red curves). In addition, if the value of $\nu_{1}$ is
increased to 2.5, then the total infected population oscillates over time
(orange curve). If the value of $\nu_{1}$ is further increased to 3.5, then
the infected population persists at a much lower level (purple curve) over
long time.

(a) (b)

Figure 10: Bifurcation diagram (endemic bubble) with respect to mutation rate
($\nu_{1}$) for the fitted values of parameters given in Table 3. The blue
color shows the upper limit of the limit cycle and red color shows the lower
limit of the cycle.

(a) (b) (c)

Figure 11: Impact of mutation rate on the infected populations $I_{1}$ and
$I_{2}$. (a) The asymptotically stable solutions $I_{1}$ and $I_{2}$ of system
(1) for the parametric values in Table 3 and $\nu_{1}=1$. (b) The periodic
solutions $I_{1}$ and $I_{2}$ of system (1) for $\nu_{1}=2.5$. (c) The
asymptotically stable solutions $I_{1}$ and $I_{2}$ of system (1) for
$\nu_{1}=3.53$. Figure 12: Variation of the total infected population with
respect to different mutation rate ($\nu_{1}$).

##  8 Sensitivity analysis

Our ultimate goal in developing a mathematical model of two strains with
imperfect vaccine and mutation is to find the role of different parameters to
control the disease. From the viewpoints of biological significance, $R_{0}$
plays a vital role in determining the severity (burden of disease), outcome
and process of the infection. This section investigates how a percentage
change in key parameters in the model affects (change) the basic reproduction
number. In the system (1), the most critical parameter to reduce is the basic
reproduction numbers of different strains. If the basic reproduction number is
brought below one, the disease infection will be eliminated. Even if the basic
reproduction number cannot be brought below one, sensitivity analysis may help
to determine which parameter, if acted upon, will bring the largest reduction
in the basic reproduction numbers. Sensitivity indices measure the percentage
change of a key quantity, such as the basic reproduction number, in response
to a percentage change of a parameter in that quantity. Sensitivity analysis
is carried out on each parameter, which is utilized to recognize and check
parameters responsible for impacting the basic reproductive number. The
normalized sensitivity indices, also called elasticity of the quantity $Q$
with respect to the parameter $p$, is defined as follows [69]:

$\epsilon_{p}^{Q}=\frac{\partial Q}{\partial p}\frac{p}{Q}.$ (23)

Elasticities can be positive or negative. A positive sign says that quantity
$Q$ increases with the increase in parameter $p$, while a negative sign says
that quantity $Q$ decreases with the increase of parameter $p$.

Table 5: List of elasticities of $R_{1}$ and $R_{2}$. Parameters | Values of elasticities of $R_{1}$ | Values of elasticities of $R_{2}$
---|---|---
$\Lambda$ | 1 | 1
$\beta_{1}$ | 1 | -
$\beta_{2}$ | - | 1
$a_{1}$ | 0.00029 | -
$a_{2}$ | - | 0.00029
$\gamma$ | 0.0014 | 0.00075
$\mu$ | -1 | -1
$p$ | -0.0014 | -0.00077
$\delta_{1}$ | 0.00048 | -
$\delta_{2}$ | - | 0.00116
$d_{1}$ | -0.0143 | -
$d_{2}$ | - | -0.00015
$\alpha_{1}$ | -0.9234 | -
$\alpha_{2}$ | - | -0.998
$\nu_{1}$ | -0.0615 | -

We compute the elasticity indices for the basic reproduction numbers $R_{1}$
and $R_{2}$ given in Table 5. In Table 5, the parameters with positive
sensitivity indices are those parameters that have a great influence on the
development of the disease in the community if their values are increasing.
This is because the basic reproduction number increases as their value
increases; that is, the average number of secondary cases of infection
increases in the community. Also, all the parameters in which their indices
are negative can curtail the infection in the community as their values
increase while the others are left constant. As their values increase, the
basic reproduction number decreases, which reduces the endemicity of the
disease in the community.

It is evident from the values of elasticities that the reproduction number
will experience the highest impact with change to recruitment rate
$(\Lambda),$ transmission rates ($(\beta_{1})$ and $(\beta_{2})$), and natural
death rate ($\mu$). In general, acting on any parameters produces a similar
change in the reproduction number. For the basic reproduction number $R_{0}$,
the elasticities concerning the parameters $\Lambda,\beta_{1},$ and
$\beta_{2},$ are fairly similar, which means these parameters have the same
impact on $R_{0}$.

We also perform sensitivity analysis using the methodology of Latin Hypercube
Sampling (LHS) and partial rank correlation coefficients (PRCCs) [49] to
investigate the dependence of $R_{0}$ on the different parameters. From Figure
13, we observe that recruitment rate ($\Lambda$), transmission rates
($\beta_{1},\beta_{2}$), natural death rate ($\mu$), and recovery rates
($\alpha_{1}$, $\alpha_{2}$) are the most sensitive parameters for $R_{0}.$ To
generate the LHS matrices, we assume that all the model parameters are
uniformly distributed. Then using the baseline values from Tables 2 and 3, a
total of 1000 simulations per LHS run are carried out.

We also examine the impact of sensitivities of the parameters on the
population size of cumulative cases $(C)$ and cumulative deaths $(D)$. From
Figure 14, we observe that $\beta_{1},a_{1},\beta_{2},a_{2}$ are the most
sensitive parameters to the cumulative cases, which means that the value of
these parameters increases the cumulative cases will increase. This result
implies that we should control these parameters to reduce the cumulative
cases. From Figure 15, we observe that $\beta_{1},a_{1},d_{1},\alpha_{1}$ are
the most sensitive parameters to the cumulative deaths and
$\beta_{1},a_{1},d_{1}$ has positive sensitivity and $\alpha_{1}$ has negative
sensitivity, indicating that we should control the parameters
$\beta_{1},a_{1},d_{1}$ and promote $\alpha_{1}$ to reduce the cumulative
deaths.

(a) (b)

Figure 13: PRCC sensitivity on $R_{0}$. Figure 14: PRCC sensitivity on
cumulative cases (C). Figure 15: PRCC sensitivity on cumulative deaths (D).

##  9 Discussion

The novel coronavirus has rapidly emerged as a disease COVID-19 and evolved as
pandemics worldwide. The emergence of new variants of SARS-CoV-2 could
complicate mitigation efforts. Reducing the transmission of the novel
coronavirus pandemic has been the massive responsibility of the intellects of
every public health agencies, Govt. officials, and millions of populations
worldwide. According to the World Health Organization (WHO) [71], over 26
crore populations are infected with the COVID-19 as of December 4, 2021. Now,
the Govt. of different countries are trying to give safeguards to the
populations via vaccination. Since every vaccine is imperfect to the disease
and the virus mutates over time, mathematical models may help understand the
dynamics of transmission and control of the novel coronavirus.

Taking care of the pandemic scenario, we proposed and analyzed an SEIR type
multi-strain mathematical model with imperfect vaccine and mutation (1). We
also assumed that recovered individuals of native strain could become infected
by the mutant strain, but they have some immunity against the mutant strain,
and it has been modeled by multiplying a reduction coefficient to the
transmission rate. The model we developed takes the form of a deterministic
system of nonlinear differential equations. The model is rigorously analyzed
to gain insight into its dynamical features. Theoretical and numerical
analysis of the proposed system (1) have been carried out using stability
theory. Positivity and boundedness of the solutions of system (1) have been
studied, and the system (1) is well-posed. We found that the DFE $(D^{0})$ is
globally stable when $R_{0}<1$ in the presence of an imperfect vaccine, i.e.,
the native and mutant strains will be eliminated in the community whenever
$R_{0}<1.$ In other words, the imperfect vaccine against COVID-19 can lead to
the eradication of the pandemic if it can bring (or maintain) $R_{0}$ to a
value less than unity. Further, we analyzed the existence and stability of the
mutant dominant equilibrium ($D^{2}$) by constructing a suitable Lyapunov
function. Global stability of coexistence (positive) equilibrium of both
native and mutant strains has been investigated for the case $\nu_{1}=0$ and
$\delta_{3}=0$, by constructing a suitable Lyapunov function. This result
implies that if the mutation rate is zero and recovered individuals from
native strain have 100% immunity against mutant strain, the disease persists
in the environment when $R_{0}>1$. Moreover, we investigated that system (1)
undergoes a Hopf bifurcation. The transcritical bifurcation was also
investigated for system (1) by using center manifold theory. By numerical
simulation, we discovered that system (1) exhibits an interesting dynamics
called an endemic bubble, which means the system losses its stability when the
mutation rate ($\nu_{1}$) is larger, and oscillations occur, and further
larger value of $\nu_{1}$, the system regains its stability (refer the Figure
4). The occurrence of Hopf bifurcation was also ensured by varying the
parameter $\delta_{3}$ (Figure 5), which implies that the disease will appear
repeatedly for the higher value of $\delta_{3}$. We also computed an
expression for herd immunity and a threshold value for the vaccination rate
that suggests disease control implications. The mentioned necessary herd
immunity (induced by vaccination) percentage may not be realistically
achievable. Indeed, Dr. Anthony Fauci, a member of the US Presidential Task
Force on COVID-19, said on June 29, 2020, that an imaginary anti-COVID-19
vaccine might not attain the requisite high immunity in the US if many people
deny getting it [11]. One way to get around this requirement for high vaccine
coverage to attain high immunity is to merge the vaccination program with
other anti-COVID-19 intervention strategies, such as social distancing, the
use of face masks in public, etc.

Furthermore, we parameterized the proposed mathematical model using the data
of the COVID-19 pandemic in India, for assessing the potential community-wide
impact of an imperfect vaccine against COVID-19 and mutation of the virus.
Using the data of cumulative confirmed cases and cumulative deaths of COVID-19
in India from March 1, 2021 to September 27, 2021. By employing the MCMC
method to fit the model (1) with the cumulative cases and cumulative deaths,
the mean values, the standard deviation, and Geweke values of the unknown
parameters were estimated.

We carried out numerical simulations to measure the population-level impact of
the vaccine’s efficacies. The results exhibited that the COVID-19 burden (as
measured in terms of cumulative cases and cumulative deaths) decreases with
increasing vaccine efficacies, as expected. We investigated how the cumulative
number of cases varies with different values of vaccine efficacies to
different strains, mutation rate, and immunity against the mutant strain in
India. Our analysis showed what would happen if the vaccine efficacy to native
strain was increased from its baseline value, as shown in Figure 7(a). The
results showed that if the vaccine efficacy $(1-\delta_{1})=90\%$ was
increased by its baseline value $(1-\delta_{1})=75\%$, then the cumulative
cases would be decreased by 3.27% by February 7, 2022. Figure 7(b)
demonstrated that if the vaccine efficacy to mutant strain increased as
$(1-\delta_{2})=75\%$ from its baseline value $(1-\delta_{2})=40\%$, then
cumulative cases would be decreased by 0.7%. We also assessed the impact of
immunity against mutant strain $(1-\delta_{3})$. We analyzed that if the value
of $(1-\delta_{3})$ is considered as 25%, which is decreased from its baseline
value $(1-\delta_{3})=90\%$, then the cumulative cases would be increased by
19.3% by February 7, 2022, shown in 7(c). We also analyzed the impact of
mutation rate on cumulative cases, which represents the increase in cumulative
cases with increasing mutation rate, showing in 7(d). Our analysis also showed
the combined impact of vaccine efficacies to both strains, which were
increased to 20%, 40%, 60%, 80%, respectively, as shown in Figure 8(a). The
result demonstrated that increasing the vaccine efficacy to the native strain
is more influential in reducing the cumulative number of cases than the
efficacy to mutant strain. We have shown in the contour plot (Figure 8(b))
that there would be higher cumulative cases for the lower immunity against
mutant strain even for the lower value of the mutation rate. There will be
fewer cumulative cases if immunity against mutant strain is high. The contour
plot (Figure 8(c)) exhibited that cumulative cases would be high for the
higher mutation rate and lower vaccine efficacy to the mutant strain. This
implies that there is a need to increase in the vaccine efficacy
$(1-\delta_{2})$ to decrease the cumulative cases of COVID-19 in the presence
of the mutation. The cumulative number of deaths has also been varied with
different values of vaccine’s efficacies to different strains in Figure 9. The
results in Figure 9(a) showed that if the vaccine efficacy
$(1-\delta_{1})=90\%$ was increased by its baseline value
$(1-\delta_{1})=75\%$, then the cumulative deaths would be decreased by 5.2%
by February 7, 2022. Figure 9(b) demonstrated that if the vaccine efficacy to
mutant strain increased as $(1-\delta_{2})=75\%$ from its baseline value
$(1-\delta_{2})=40\%$, then cumulative deaths would be decreased by 1.1%.

The endemic bubble phenomenon has been verified for COVID-19 by plotting the
bifurcation diagram and periodic solutions for the fitted values from Table 3
in Figures 10 and 11, respectively. These results imply that for a range of
values of mutation rate ($\nu_{1}$), the COVID-19 outbreak will occur
repeatedly. However, for the lower and higher enough values of $\nu_{1}$, the
disease remains stable in the environment. These dynamics represent a worse
scenario for the mutation of the virus in the sense of the occurrence of
periodic solutions. Furthermore, the total infected population has also been
plotted for different values of mutation rate $\nu_{1}$ and a wondering
dynamics has been observed. This result shows that for a high rate of
mutation, the total infected population would be stable at a lower level over
a long time. However, oscillation occurs for an intermediate range of mutation
rate, implying that mutation may be worse in the sense of oscillation but
better in the sense of stability at a lower infection level over a long time.

Furthermore, sensitivity analysis was performed to reveal the relative
significance of the key epidemiological parameters of the system (1), which
are $\Lambda,\beta_{1},\beta_{2},\alpha_{1},\alpha_{2}$, because these
parameters should be given priority to effectively control the disease. The
PRCCs of the reproduction number (Figure 13) showed that reducing
$\Lambda,\beta_{1},\beta_{2}$, which may be realized by strong control
measures, such as lockdown, using face mask, travel restriction, isolation of
infected individuals, contact tracing, can significantly reduce $R_{0}$ and
thus lower the transmission risk of COVID-19. The PRCCs of the cumulative
cases and deaths suggest that the parameters $\beta_{1},a_{1},\beta_{2},$ and
$a_{2}$ should be controlled in order to reduce the cumulative cases, and
$\beta_{1},a_{1},$ and $d_{1}$ should be controlled in order to reduce the
cumulative deaths.

We have to point out that the results in contour plots revealed an interesting
problem, i.e., there may exist an optimal strategy of supplementing the
vaccine with different efficacies, in order to ensure that the fewer
cumulative cases and cumulative deaths, and that the expenditure is
economical, a topic for our future work. For preserving public health, it must
be focused on how the facility of medical resources availability affects the
pandemic of COVID-19, which could be included in our model as future work.
Formulating and analyzing the model with other non-pharmaceutical
interventions [73] with vaccination will be more helpful to control the spread
of multi-strain disease. For instance, Iboi et al. [33] proposed a dynamical
model to assess the impact of vaccine strategy with other public health
intervention strategies; in particular, the authors assessed the impact of
vaccine efficacy combined with mask efficacy concerning COVID-19 in the US.
Their results confirmed that the elimination of COVID-19 is more feasible if
the vaccine program is combined with other interventions. We focused on the
situation in India, but the model can be extended to describe the efficacies
of vaccines to different strains with mutation rates and developed immunity
against mutant strains in other countries.

## References

  * [1] Ahmed SE (2008) Markov chain monte carlo: Stochastic simulation for bayesian inference. https://doi.org/10.1198/tech.2008.s542.
  * [2] Alexander ME, Bowman C, Moghadas SM, Summers R, Gumel AB, Sahai BM (2004) A vaccination model for transmission dynamics of influenza. SIAM Journal on Applied Dynamical Systems 3(4):503-524. https://doi.org/10.1137/030600370.
  * [3] Anderson RM (1992) The concept of herd immunity and the design of community-based immunization programmes. Vaccine 10:928e935. https://doi.org/10.1016/0264-410X(92)90327-G.
  * [4] Anderson RM, May RM (1985) Vaccination and herd immunity to infectious diseases. Nature 318:323e329. https://doi.org/10.1038/318323a0.
  * [5] Arino J, McCluskey CC, van den Driessche P (2003) Global results for an epidemic model with vaccination that exhibits backward bifurcation. SIAM Journal on Applied Mathematics 64(1):260-276. https://doi.org/10.1137/S0036139902413829.
  * [6] Arruda EF, Pastore DH, Dias CM, Das SS (2021) Modelling and Optimal Control of Multi Strain Epidemics, with Application to COVID-19. arXiv preprint 2101:08137. 10.1371/journal.pone.0257512.
  * [7] Bajiya VP, Bugalia S, Tripathi JP (2020) Mathematical modeling of COVID-19: impact of non-pharmaceutical interventions in India. Chaos: An Interdisciplinary Journal of Nonlinear Science 30(11):113143. https://doi.org/10.1063/5.0021353.
  * [8] Bugalia S, Bajiya VP, Tripathi JP, Li MT, Sun GQ (2020) Mathematical modeling of COVID-19 transmission: the roles of intervention strategies and lockdown. Mathematical Biosciences and Engineering 17(5):5961-5986. 10.3934/mbe.2020318.
  * [9] Bugalia S, Tripathi JP, Wang H (2021) Mathematical modeling of intervention and low medical resource availability with delays: Applications to COVID-19 outbreaks in Spain and Italy. Mathematical Biosciences and Engineering 18(5):5865-5920. 10.3934/mbe.2021295.
  * [10] Bussiness Insider (2021) https://tinyurl.com/6q3g47b8. Accessed 2 November 2021.
  * [11] Cable News Network (2020) Fauci says covid-19 vaccine may not get us to herd immunity if too many people refuse to get it, CNN. https://edition.cnn.com/2020/06/28/health/fauci-coronavirus-vaccine-contact-tracing-aspen/index.html. Accessed 10 September 2021.
  * [12] Cai L, Xiang J, Li X, Lashari AA (2012) A two-strain epidemic model with mutant strain and vaccination. Journal of Applied Mathematics and Computing 40(1):125-142. https://doi.org/10.1007/s12190-012-0580-x.
  * [13] Cai LM, Li Z, Song X (2018) Global analysis of an epidemic model with vaccination. Journal of Applied Mathematics and Computing 57(1):605-628. https://doi.org/10.1007/s12190-017-1124-1.
  * [14] Castillo-Chavez C, Blower S, Van den Driessche P, Kirschner D, Yakubu AA, editors (2002) Mathematical approaches for emerging and reemerging infectious diseases: models, methods, and theory. Springer Science & Business Media.
  * [15] Castillo-Chavez C, Feng Z, Huang W (2002) On the computation of $R_{0}$ and its role in global stability. IMA Volumes in Mathematics and Its Applications 125:229-250.
  * [16] Castillo-Chavez C, Song B (2004) Dynamical models of tuberculosis and their applications. Mathematical Biosciences & Engineering 1(2):361. 10.3934/mbe.2004.1.361.
  * [17] Centers for Disease Control and Prevention (CDC) (2021). https://www.cdc.gov/coronavirus/2019-ncov/variants/variant.html. Accessed 13 October 2021.
  * [18] Centers for Disease Control and Prevention (CDC) (2021). https://www.cdc.gov/coronavirus/2019-ncov/your-health/reinfection.html. Accessed 13 October 2021.
  * [19] Deng X, Garcia-Knight MA, Khalid MM, Servellita V, Wang C, Morris MK, Sotomayor-González A, Glasner DR, Reyes KR, Gliwa AS, Reddy NP (2021) Transmission, infectivity, and neutralization of a spike L452R SARS-CoV-2 variant. Cell. https://doi.org/10.1016/j.cell.2021.04.025.
  * [20] Diekmann O, Heesterbeek JA, Roberts MG (2010) The construction of next-generation matrices for compartmental epidemic models. Journal of the Royal Society Interface 7(47):873-885. https://doi.org/10.1098/rsif.2009.0386.
  * [21] Elbasha EH, Gumel AB (2021) Vaccination and herd immunity thresholds in heterogeneous populations. Journal of Mathematical Biology 83(6-7). https://doi.org/10.1007/s00285-021-01686-z.
  * [22] Eron JJ, Vernazza PL, Johnston DM, Seillier-Moiseiwitsch F, Alcorn TM, Fiscus SA, Cohen MS (1998) Resistance of HIV-1 to antiretroviral agents in blood and seminal plasma: implications for transmission. Aids 12(15):F181-F189.
  * [23] Friedman TL, Is Sweden doing it right? New York Times, https://www.nytimes.com/2020/04/28/opinion/coronavirus-sweden.html. Accessed 3 October 2021.
  * [24] Fudolig M, Howard R (2020) The local stability of a modified multi-strain SIR model for emerging viral strains. PloS One 15(12):e0243408. https://doi.org/10.1371/journal.pone.0243408.
  * [25] Gonzalez-Parra G, Martínez-Rodríguez D, Villanueva-Micó RJ (2021) Impact of a new SARS-CoV-2 variant on the population: A mathematical modeling approach. Mathematical and Computational Applications 26(2):25. https://doi.org/10.3390/mca26020025.
  * [26] Gukenheimer J, Holmes P (1983) Nonlinear oscillations, Dynamical Systems, and Bifurcation of Vector Fields. Springer-Verlag NY.
  * [27] Gumel AB, McCluskey CC, Watmough J (2006) An SVEIR model for assessing potential impact of an imperfect anti-SARS vaccine. Mathematical Biosciences & Engineering 3(3):485. 10.3934/mbe.2006.3.485.
  * [28] Gupta S, Ferguson NM, Anderson RM (1997) Vaccination and the population structure of antigenically diverse pathogens that exchange genetic material. Proceedings of the Royal Society of London. Series B: Biological Sciences 264(1387):1435-1443. https://doi.org/10.1098/rspb.1997.0200.
  * [29] Haario H, Laine M, Mira A, Saksman E (2006) DRAM: efficient adaptive MCMC. Statistics and computing 16(4):339-354. https://doi.org/10.1007/s11222-006-9438-0.
  * [30] Haario H, Saksman E, Tamminen J (2001) An adaptive Metropolis algorithm. Bernoulli 7(2):223-242. https://doi.org/10.2307/3318737.
  * [31] Health, The Sciences, https://science.thewire.in/the-sciences/covid-19-reinfection-hong-kong-man-immunity/. Accessed 10 November 2021.
  * [32] Hethcote HW (2000) The mathematics of infectious diseases. SIAM review Society for Industrial and Applied Mathematics 42(4):599-653. https://doi.org/10.1137/S0036144500371907.
  * [33] Iboi EA, Ngonghala CN, Gumel AB (2020) Will an imperfect vaccine curtail the COVID-19 pandemic in the US? Infectious Disease Modelling 5:510-524. https://doi.org/10.1016/j.idm.2020.07.006.
  * [34] Iboi EA, Sharomi OO, Ngonghala CN, Gumel AB (2020) Mathematical modeling and analysis of COVID-19 pandemic in Nigeria. Mathematical Biosciences and Engineering 17(6):7192-7220. 10.3934/mbe.2020369.
  * [35] Iwami S, Takeuchi Y, Liu X (2007) Avian–human influenza epidemic model. Mathematical Biosciences 207(1):1-25. https://doi.org/10.1016/j.mbs.2006.08.001.
  * [36] Johns Hopkins Medicine, https://www.hopkinsmedicine.org/health/conditions-and-diseases/coronavirus/a-new-strain-of-coronavirus-what-you
-should-know. Accessed 17 September 2021.
  * [37] Kermack WO, McKendrick AG (1927) A contribution to the mathematical theory of epidemics. Proceedings of the Royal Society of London, Series A, Containing papers of a mathematical and physical character 115(772):700-721. https://doi.org/10.1098/rspa.1927.0118.
  * [38] Khyar O, Allali K (2020) Global dynamics of a multi-strain SEIR epidemic model with general incidence rates: application to COVID-19 pandemic. Nonlinear Dynamics 102(1):489-509. https://doi.org/10.1007/s11071-020-05929-4.
  * [39] Korber B, Fischer WM, Gnanakaran S, Yoon H, Theiler J, Abfalterer W, Hengartner N, Giorgi EE, Bhattacharya T, Foley B, Hastie KM (2020) Tracking changes in SARS-CoV-2 Spike: evidence that D614G increases infectivity of the COVID-19 virus. Cell 182(4):812-827. https://doi.org/10.1016/j.cell.2020.06.043.
  * [40] LaSalle JP (1976) The stability of dynamical systems. Society for Industrial and Applied Mathematics.
  * [41] Lauer SA, Grantz KH, Bi Q, Jones FK, Zheng Q, Meredith HR, Azman AS, Reich NG, Lessler J (2020) The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: estimation and application. Annals of internal medicine 172(9):577-582. https://doi.org/10.7326/M20-0504.
  * [42] Lemieux JE, Li JZ (2021) Uncovering Ways that Emerging SARS-CoV-2 Lineages May Increase Transmissibility. The Journal of Infectious Diseases 223(10):1663–1665. https://doi.org/10.1093/infdis/jiab083.
  * [43] Li J, Zhou Y, Ma Z, Hyman JM (2004) Epidemiological models for mutating pathogens. SIAM Journal on Applied Mathematics 65(1):1-23. https://doi.org/10.1137/S0036139903430185.
  * [44] Liu L, Ren X, Liu X (2018) Dynamical behaviors of an influenza epidemic model with virus mutation. Journal of Biological Systems 26(03):455-472. https://doi.org/10.1142/S0218339018500201.
  * [45] Liu M, Liz E, Rost G (2015) Endemic bubbles generated by delayed behavioral response: global stability and bifurcation switches in an SIS model. SIAM Journal on Applied Mathematics 75(1):75-91. https://doi.org/10.1137/140972652.
  * [46] Liu WM (1994) Criterion of Hopf bifurcations without using eigenvalues. Journal of Mathematical Analysis and Applications 182(1):250-256. https://doi.org/10.1006/jmaa.1994.1079.
  * [47] Liu X, Takeuchi Y, Iwami S (2008) SVIR epidemic models with vaccination strategies. Journal of Theoretical Biology 253(1):1-11. https://doi.org/10.1016/j.jtbi.2007.10.014.
  * [48] MacIntyre CR, Costantino V, Trent MJ (2021) Modelling of COVID-19 vaccination strategies and herd immunity, in scenarios of limited and full vaccine supply in NSW, Australia. Vaccine. https://doi.org/10.1016/j.vaccine.2021.04.042.
  * [49] Marino S, Hogue IB, Ray CJ, Kirschner DE (2008) A methodology for performing global uncertainty and sensitivity analysis in systems biology. Journal of Theoretical Biology 254(1):178-196. https://doi.org/10.1016/j.jtbi.2008.04.011.
  * [50] Martcheva M (2015) An introduction to mathematical epidemiology. New York: Springer. https://doi.org/10.1007/978-1-4899-7612-3.
  * [51] Martcheva M, Iannelli M, Li XZ (2007) Subthreshold coexistence of strains: the impact of vaccination and mutation. Mathematical Biosciences and Engineering 4(2):287. 10.3934/mbe.2007.4.287.
  * [52] May RM, Nowak MA (1995) Coinfection and the evolution of parasite virulence. Proceedings of the Royal Society of London. Series B, Biological Sciences 261(1361):209-215. https://doi.org/10.1098/rspb.1995.0138.
  * [53] May RM, Nowak MA (1994) Superinfection, metapopulation dynamics, and the evolution of diversity. Journal of Theoretical Biology 170(1):95-114. https://doi.org/10.1006/jtbi.1994.1171.
  * [54] McLean AR (1995) Vaccination, evolution and changes in the efficacy of vaccines: a theoretical framework. Proceedings of the Royal Society of London. Series B, Biological Sciences 261(1362):389-393. https://doi.org/10.1098/rspb.1995.0164.
  * [55] Ministry of Health and Family Welfare, Government of India. https://www.mohfw.gov.in/covid_vaccination/vaccination/faqs.html#what-to-expect-before-vaccination-2. Accessed 16 October 2021.
  * [56] Nature news, https://www.nature.com/articles/d41586-021-01059-y. Accessed 18 October 2021.
  * [57] Olliaro P, Torreele E, Vaillant M (2021) COVID-19 vaccine efficacy and effectiveness—the elephant (not) in the room. The Lancet Microbe. https://doi.org/10.1016/S2666-5247(21)00069-0.
  * [58] Palese P, Young JF (1982) Variation of influenza A, B, and C viruses. Science 215(4539):1468-1474. 10.1126/science.7038875.
  * [59] Parton R, Hall E, Wardlaw AC (1994) Responses to Bordetella pertussis mutant strains and to vaccination in the coughing rat model of pertussis. Journal of Medical Microbiology 40(5):307-312. https://doi.org/10.1099/00222615-40-5-307.
  * [60] Porco TC, Blower SM (1998) Designing HIV vaccination policies: subtypes and cross-immunity. Interfaces 28(3):167-190. https://doi.org/10.1287/inte.28.3.167.
  * [61] Porco TC, Blower SM (2000) HIV vaccines: the effect of the mode of action on the coexistence of HIV subtypes. Mathematical Population Studies 8(2):205-229. https://doi.org/10.1080/08898480009525481.
  * [62] Public Health England (2021) investigation-of-novel-sars-cov-2-variant-variant-of-concern-20201201. https://www.gov.uk/government/publications/investigation-of-novel-sars-cov-2-variant-variant-of-concern-20201201. Accessed 27 September 2021.
  * [63] Rahimi F, Abadi AT (2021) Implications of the Emergence of a New Variant of SARS-CoV-2, VUI-202012/01. Archives of Medical Research 52(5):569–571. https://doi.org/10.1016/j.arcmed.2021.01.001.
  * [64] Sansonetti PJ, Arondel J (1989) Construction and evaluation of a double mutant of Shigella flexneri as a candidate for oral vaccination against shigellosis. Vaccine 7(5):443-450. https://doi.org/10.1016/0264-410X(89)90160-6.
  * [65] Sato S, Suzuki K, Akahane Y, Akamatsu K, Akiyama K, Yunomura K, Tsuda F, Tanaka T, Okamoto H, Miyakawa Y, Mayumi M (1995) Hepatitis B virus strains with mutations in the core promoter in patients with fulminant hepatitis. Annals of Internal Medicine 122(4):241-248. https://doi.org/10.7326/0003-4819-122-4-199502150-00001.
  * [66] Scherer A, McLean A (2002) Mathematical models of vaccination. British Medical Bulletin 62(1):187-199. https://doi.org/10.1093/bmb/62.1.187.
  * [67] Shastri J, Parikh S, Aggarwal V, Agrawal S, Chatterjee N, Shah R, Devi P, Mehta P, Pandey R (2021) Severe SARS-CoV-2 breakthrough reinfection with Delta variant after recovery from breakthrough infection by Alpha variant in a fully vaccinated health worker. Frontiers in Medicine 1379. https://doi.org/10.3389/fmed.2021.737007.
  * [68] Van den Driessche P, Watmough J (2002) Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical Biosciences 180(1-2):29-48. https://doi.org/10.1016/S0025-5564(02)00108-6.
  * [69] Martcheva M (2014) Avian flu: modeling and implications for control. Journal of Biological Systems 22(01):151-75. https://doi.org/10.1142/S0218339014500090.
  * [70] World Health Organization, https://covid19.who.int/region/searo/country/in. Accessed 29 October 2021.
  * [71] World Health Organization, https://covid19.who.int/. Accessed 29 October 2021.
  * [72] World Health Organization, WHO advisory committee on variola virus research: report of the thirteenth meeting report (World Health Organization, 2011), https://apps.who.int/iris/handle/10665/70778. Accessed 29 October 2021.
  * [73] Xue L, Jing SL, Wang H (2022) Evaluating the impacts of non-pharmaceutical interventions on the transmission dynamics of COVID-19 in Canada based on mobile network. PloS One (in press).

## Appendix A

$\displaystyle l_{1}=$
$\displaystyle-(a_{11}+b_{12}+c_{13}+d_{13}+e_{12}+f_{13}+g_{13}),$
$\displaystyle l_{2}=$ $\displaystyle-
a_{12}b_{11}+c_{13}d_{13}+c_{13}e_{12}-c_{14}e_{11}+d_{13}e_{12}-d_{14}f_{11}+(c_{13}+d_{13}+e_{12})f_{13}$
$\displaystyle+g_{13}(c_{13}+d_{13}+e_{12}+f_{13})+b_{12}(c_{13}+d_{13}+e_{12}+f_{13}+g_{13})$
$\displaystyle+a_{11}(b_{12}+c_{13}+d_{13}+e_{12}+f_{13}+g_{13}),$
$\displaystyle l_{3}=$ $\displaystyle
a_{11}(-b_{12}(c_{13}+d_{13}+e_{12}+f_{13}+g_{13})-c_{13}(d_{13}+e_{12}+f_{13}+g_{13})+c_{14}e_{11}$
$\displaystyle-
g_{13}(d_{13}+e_{12}+f_{13})-d_{13}e_{12}-d_{13}f_{13}+d_{14}f_{11}-e_{12}f_{13})$
$\displaystyle+a_{12}b_{11}(c_{13}+d_{13}+e_{12}+f_{13}+g_{13})-g_{13}(b_{12}(c_{13}+d_{13}+e_{12}+f_{13})$
$\displaystyle+c_{13}(d_{13}+e_{12}+f_{13})-c_{14}e_{11}+f_{13}(d_{13}+e_{12})+d_{13}e_{12}-d_{14}f_{11})$
$\displaystyle-
e_{11}(a_{13}c_{11}+b_{13}c_{12})-c_{13}d_{13}e_{12}+c_{13}d_{14}f_{11}-c_{13}d_{13}f_{13}-c_{13}e_{12}f_{13}$
$\displaystyle-
f_{11}(a_{14}d_{11}+b_{14}d_{12})+d_{14}e_{12}f_{11}-d_{13}e_{12}f_{13}+c_{14}e_{11}(d_{13}+f_{13})-d_{15}f_{11}g_{12}$
$\displaystyle-
b_{12}(c_{13}(d_{13}+e_{12}+f_{13})-c_{14}e_{11}+f_{13}(d_{13}+e_{12})+d_{13}e_{12}-d_{14}f_{11}),$
$\displaystyle l_{4}=$ $\displaystyle
a_{13}e_{11}(c_{11}(b_{12}+d_{13})-b_{11}c_{12})+d_{13}(b_{12}c_{13}e_{12}-b_{12}c_{14}e_{11}+b_{13}c_{12}e_{11})$
$\displaystyle+f_{11}(a_{14}(d_{11}(b_{12}+c_{13}+e_{12})-b_{11}d_{12})-d_{14}(e_{12}(b_{12}+c_{13})+b_{12}c_{13}-c_{14}e_{11})$
$\displaystyle+b_{14}d_{12}(c_{13}+e_{12}))-a_{12}(b_{11}(c_{13}(d_{13}+e_{12})-c_{14}e_{11}+d_{13}e_{12}-d_{14}f_{11})+b_{13}c_{11}e_{11}$
$\displaystyle+b_{14}d_{11}f_{11})-e_{11}f_{12}(a_{14}c_{11}+b_{14}c_{12})+b_{12}c_{13}d_{13}f_{13}+f_{13}(a_{13}c_{11}e_{11}$
$\displaystyle+e_{12}(d_{13}(b_{12}+c_{13})+b_{12}c_{13})-c_{14}e_{11}(b_{12}+d_{13})+b_{13}c_{12}e_{11})$
$\displaystyle-
a_{12}b_{11}f_{13}(c_{13}+d_{13}+e_{12})+d_{15}f_{11}g_{12}(b_{12}+c_{13}+e_{12})-a_{12}b_{11}g_{13}(c_{13}$
$\displaystyle+d_{13}+e_{12}+f_{13})+g_{13}(a_{13}c_{11}e_{11}+a_{14}d_{11}f_{11}+b_{12}(c_{13}(d_{13}+e_{12}+f_{13})$
$\displaystyle-
c_{14}e_{11}+f_{13}(d_{13}+e_{12})+d_{13}e_{12}-d_{14}f_{11})+b_{13}c_{12}e_{11}+b_{14}d_{12}f_{11}$
$\displaystyle+c_{13}d_{13}e_{12}+c_{13}d_{13}f_{13}-c_{13}d_{14}f_{11}+c_{13}e_{12}f_{13}-c_{14}e_{11}(d_{13}+f_{13})$
$\displaystyle+d_{13}e_{12}f_{13}-d_{14}e_{12}f_{11})+a_{11}(b_{13}c_{12}e_{11}+b_{14}d_{12}f_{11}+c_{13}(f_{13}(d_{13}+e_{12})$
$\displaystyle+d_{13}e_{12}-d_{14}f_{11})+d_{13}e_{12}f_{13}-d_{14}e_{12}f_{11}+d_{15}f_{11}g_{12})+a_{11}(b_{12}(c_{13}(d_{13}$
$\displaystyle+e_{12}+f_{13}+g_{13})-c_{14}e_{11}+d_{13}(e_{12}+f_{13}+g_{13})-d_{14}f_{11}+g_{13}(e_{12}+f_{13})$
$\displaystyle+e_{12}f_{13})+g_{13}(c_{13}(d_{13}+e_{12}+f_{13})+d_{13}(e_{12}+f_{13})-d_{14}f_{11}+e_{12}f_{13})$
$\displaystyle-c_{14}e_{11}(d_{13}+f_{13}+g_{13})),$ $\displaystyle l_{5}=$
$\displaystyle
a_{12}(b_{11}(g_{13}(c_{13}(d_{13}+e_{12}+f_{13})+f_{13}(d_{13}+e_{12})+d_{13}e_{12}-d_{14}f_{11})+c_{13}d_{13}e_{12}$
$\displaystyle+c_{13}d_{13}f_{13}-c_{13}d_{14}f_{11}+c_{13}e_{12}f_{13}-c_{14}e_{11}(d_{13}+f_{13}+g_{13})+d_{13}e_{12}f_{13}$
$\displaystyle-
d_{14}e_{12}f_{11}+d_{15}f_{11}g_{12})+b_{13}c_{11}e_{11}(d_{13}+f_{13}+g_{13})+b_{14}(d_{11}f_{11}(c_{13}+e_{12}+g_{13})$
$\displaystyle-
c_{11}e_{11}f_{12}))-a_{13}e_{11}(-b_{11}c_{12}(d_{13}+f_{13}+g_{13})+b_{12}c_{11}(d_{13}+f_{13}+g_{13})$
$\displaystyle+c_{11}(g_{13}(d_{13}+f_{13})+d_{13}f_{13}-d_{14}f_{11}))-g_{13}(f_{11}(a_{14}(d_{11}(c_{13}+e_{12})-b_{11}d_{12})$
$\displaystyle+b_{14}d_{12}(c_{13}+e_{12})-c_{13}d_{14}e_{12}+c_{14}d_{14}e_{11})-b_{12}(-a_{14}d_{11}f_{11}-c_{13}d_{13}e_{12}$
$\displaystyle+d_{14}f_{11}(c_{13}+e_{12})+c_{14}d_{13}e_{11})-e_{11}f_{12}(a_{14}c_{11}+b_{14}c_{12})+b_{13}c_{12}d_{13}e_{11}$
$\displaystyle+f_{13}(b_{12}(c_{13}(d_{13}+e_{12})-c_{14}e_{11}+d_{13}e_{12})+b_{13}c_{12}e_{11}+c_{13}d_{13}e_{12}-c_{14}d_{13}e_{11})$
$\displaystyle+a_{11}(b_{12}(c_{13}(d_{13}+e_{12}+f_{13})-c_{14}e_{11}+f_{13}(d_{13}+e_{12})+d_{13}e_{12}-d_{14}f_{11})$
$\displaystyle+b_{13}c_{12}e_{11}+b_{14}d_{12}f_{11}+c_{13}d_{13}e_{12}+c_{13}d_{13}f_{13}-c_{13}d_{14}f_{11}+c_{13}e_{12}f_{13}$
$\displaystyle-
c_{14}e_{11}(d_{13}+f_{13})+d_{13}e_{12}f_{13}-d_{14}e_{12}f_{11}))+a_{14}e_{11}f_{12}(c_{11}(b_{12}+d_{13})-b_{11}c_{12})$
$\displaystyle+a_{14}f_{11}(b_{11}d_{12}(c_{13}+e_{12})-b_{12}d_{11}(c_{13}+e_{12})-c_{13}d_{11}e_{12}+c_{14}d_{11}e_{11})$
$\displaystyle-
d_{15}f_{11}g_{12}(b_{12}(c_{13}+e_{12})+c_{13}e_{12}-c_{14}e_{11})+b_{14}c_{12}d_{13}e_{11}f_{12}+f_{11}(d_{14}(b_{12}c_{13}e_{12}$
$\displaystyle-
b_{12}c_{14}e_{11}+b_{13}c_{12}e_{11})+b_{14}d_{12}(c_{14}e_{11}-c_{13}e_{12}))-d_{13}f_{13}(b_{12}c_{13}e_{12}-b_{12}c_{14}e_{11}$
$\displaystyle+b_{13}c_{12}e_{11})+a_{11}(-f_{13}(b_{12}(c_{13}(d_{13}+e_{12})-c_{14}e_{11}+d_{13}e_{12})+b_{13}c_{12}e_{11}+c_{13}d_{13}e_{12}$
$\displaystyle-
c_{14}d_{13}e_{11})-b_{12}c_{13}d_{13}e_{12}+b_{12}c_{13}d_{14}f_{11}-d_{15}f_{11}g_{12}(b_{12}+c_{13}+e_{12})+b_{12}c_{14}d_{13}e_{11}$
$\displaystyle+b_{12}d_{14}e_{12}f_{11}-b_{13}c_{12}d_{13}e_{11}+b_{14}c_{12}e_{11}f_{12}-b_{14}c_{13}d_{12}f_{11}-b_{14}d_{12}e_{12}f_{11}$
$\displaystyle+c_{13}d_{14}e_{12}f_{11}-c_{14}d_{14}e_{11}f_{11}),$
$\displaystyle l_{6}=$ $\displaystyle
a_{13}e_{11}(-g_{13}(b_{11}c_{12}(d_{13}+f_{13})-c_{11}d_{13}f_{13}+c_{11}d_{14}f_{11})-b_{11}c_{12}d_{13}f_{13}+b_{11}c_{12}d_{14}f_{11}$