Qualitative and Quantitative Analyses of COVID-19 Dynamics

Abstract

COVID-19 is a highly contagious disease which has spread across the world. A deterministic model that considers an important component of individuals with vertically transmitted underlying diseases (high-risk susceptible individuals), rather than the general public, is formulated in this paper. We also consider key parameters that are concerned with the disease. An epidemiological threshold, $R_0$, is computed using next-generation matrix approach. This is used to establish the existence and global stability of equilibria. We identify the most sensitive parameters which effectively contribute to change the disease dynamics with the help of sensitivity analysis. Our results reveal that increasing contact tracing of the exposed individuals who are tested for COVID-19 and hospitalizing them, largely has a negative impact on $R_0$. Results further reveal that transmission rate between low-risk/high-risk susceptible individuals and symptomatic infectious individuals $\beta$ and incubation rate of the exposed individuals $\sigma$ have positive impact on $R_0$. Numerical simulations show that there are fewer high-risk susceptible individuals than the general public when $R_0<1$. This may be due to the fact that high-risk susceptible individuals may prove a bit more difficult to control than the low-risk susceptible individuals as a result of inherited underlying diseases present in them. We thus conclude that high level of tracing and hospitalizing the exposed individuals, as well as adherence to standard precautions and wearing appropriate Personal Protective Equipment (PPE) while handling emergency cases, are needed to flatten the epidemic curve.

1. Introduction

Until now, the COVID-19 pandemic remains a global concern due to the continuous increase in the number of infectious individuals. The number of confirmed cases has been growing very fast on a daily basis and it has been declared a worldwide pandemic disease. The disease started in Wuhan, Hubei Province, China, in 2020 and since then, it has spread across the world [1]. Presently, the pandemic has disrupted economies all over the world. Since December 2019, the world has reported 148,329,348 confirmed cases and 3,128,962 death cases [2]. High fever, severe chest pain, body aches, headache, difficulty in breathing, fatigue, pains, sore throat, skin rashes, running nose, taste loss, and diarrhea are the symptoms of COVID-19, but the main and most common symptoms are dry cough and fatigue [1].

With ongoing community transmission from asymptomatic infectious individuals, disease burden is expected to rise. Consequently, there will be an ongoing need for people with underlying diseases who are front-line health care workers in patient-facing roles. This is due to the fact that their work requires close personal exposure to patients with SARS-CoV-2. They are at high risk of infection contributing to further spread [3]. SARS-Cov-2 infects all susceptible individuals. However, evidence to date suggests that a group of people is at a higher risk of contacting the virus quickly than the general public, due to compromised immune system. These are people with underlying medical conditions. An underlying disease affects the immunity, nutritional status and general well-being of patients-factors that play a critical role in protecting individuals against COVID-19. The World Health Organization emphasizes that this group of individuals (high-risk susceptible individuals) must protect themselves from COVID-19 in order to protect others. Furthermore, the highly susceptible individuals experience severe symptoms that may lead to death when they contact COVID-19. In response to the pandemic, government authorities at various levels have put in place various control measures such as imposing strict mandatory lockdown, use of a face mask in public, regular hand washing, and social distancing. Therefore, in order to reduce COVID-19 spread, contact tracing of suspected infected individuals has been stepped up in several countries and detected cases are quickly placed in isolation for prompt treatment [4].

The use of mathematical models to describe the transmission and spread of COVID-19, so as to gain insights into the disease behaviour and develop strategies for its curtailment, has been extensively studied by researchers. (see in [5,6,7,8,9,10,11,12,13,14,15,16,17,18] and the references cited therein.) Our model has a unique difference from the models developed in [10,18] as follows: the authors of [18] developed a mathematical model of COVID-19 to investigate the impact of non-pharmaceutical interventions for possible control of the disease. In their model, they only considered symptomatically infected and asymptomatically infected individuals without considering some individuals in the susceptible population who are more vulnerable to the disease than the rest of the susceptible population. These are front-line health care workers, pregnant women, children and family of COVID-19 patients. For instance, any COVID-19 patient who visits hospital for medical treatment will be attended to by these health care workers and this puts them in higher risk of contacting the virus. Classifying these individuals as a risk and treating them as a genuine threat is a national priority. Furthermore, in [10], they modeled the effects of non-pharmaceutical interventions on COVID-19 spread in Kenya. Their model excludes the class of high risk susceptible individuals which is taken into account in our model. The work in [19] investigated a class of SEIRS epidemic models with a general nonlinear incidence function. Systems of differential equations are also used as HIV infection model with transitions between uninfected cells and infected cells. The classical HIV infection divided the cells population into following three main compartments depending on cell status, target cells, infected cells producing viruses, matured virus particles, and its concentrations are denoted by $x\left(t\right),y\left(t\right)$, and $v\left(t\right)$, respectively. The authors of [20] developed a nonlinear mathematical model of MERS-COV to study the dynamical behaviour of the disease with two discrete-time delays. We highlight the main objectives of the study to achieve our aim in what follows.

(1)

It is noteworthy to mention that the high-risk susceptible individuals that we refer to, in this paper, are individuals with vertically transmitted or inherited underlying diseases, i.e., HIV, asthma and so on. These individuals are at higher-risk of contacting COVID-19 than the rest of individuals in the population and thus the need for the inclusion of high-risk susceptible individuals in order to understand the disease dynamics.

(2)

Global stability analyses of the disease-free and endemic equilibrium points are derived

(3)

Robust sensitivity analysis of the model is performed to determine the contributory effects of the factors/parameters on the spread of COVID-19.

The rest of the paper is organized as follows. Section 2 has the model formulation and computation of basic reproduction number. In Section 3, global stability of disease-free and endemic equilibrium points are established. Sensitivity analysis of the model parameters is performed in Section 4. Some simulations are done in order to confirm our theoretical results, regarding the global stability of equilibriums, in Section 5. Section 6 wraps the modeling work with conclusion.

2. Materials and Methods

The total population size, N, is divided into eight mutually exclusive compartments; S, J, E, A, B, H, C, and R, which is interpreted as (i) low-risk susceptible individuals, (ii) high-risk susceptible individuals, (iii) exposed individuals, (iv) asymptomatic infectious individuals, (v) symptomatic infectious individuals, (vi) hospitalized individuals, (vii) death class, and (viii) recovered individuals. People with vertically transmitted or inherited disease are assumed as the high-risk susceptible individuals. Individuals in S and J can be exposed to Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) through contact with either A or B, and we assume that the disease is faster transmitted in symptomatic infectious population than the asymptomatic infectious population, and include a variability factor $\pi$. It is assumed that the transmission rate in classes S and J are the same. We also assume that individuals in H are in isolation and therefore cannot transmit the disease. We assume that new recruits enter human population by birth, immigration, or emigration at the rate $\mathrm{\Lambda }$. Furthermore, there is vertical transmission recruitment into class J at a rate $\varphi$. With contact tracing, a fraction, p, of the exposed individuals is tested for COVID-19 at a rate q and become hospitalized. A fraction $ϵ$ of the remaining fraction $(1-p)$ of the exposed individuals who show no clinical symptoms of COVID-19 at the end of incubation period become asymptomatic at a progression rate $\sigma$, while the remaining fraction $(1-ϵ)(1-p)$ of the exposed individuals who show clinical symptoms at the end of incubation period becomes symptomatic at a progression rate $\sigma$. Fraction of low-risk and high-risk susceptible individuals who adhere to COVID-19 prevention guidelines (wearing of face masks, hand washing, and maintaining social distancing) is denoted as $0<\psi \le 1$, while $0<g\le 1$ represents the expected decrease in the risk of infection due to strict and correct compliance to COVID-19 prevention guidelines by low-risk and high-risk susceptible individuals. Further, a fraction of low-risk and high-risk susceptible individuals who are on lockdown within a population is denoted as $0<\nu \le 1$, while $0<m\le 1$ is the expected decrease in the risk of infection due to lock-down. The asymptomatic and symptomatic infectious individuals are hospitalized at the rates $\kappa$ and k and recover at the rates $\tau$ and d. The hospitalized individuals recover at the rate $\phi$ but die due to COVID-19 at the rate u. The symptomatic infectious individuals also die due to COVID-19 at the rate $\delta$. Dead individuals are buried at the rate a. $\eta$ and $\omega$ are the transmission rate of J and A and J and B while $\alpha$ and $\omega$ are the transmission rate of S and A and S and B. We ignore the transmission route from death compartment as the likelihood of transmission when handling human remains is low [21].

Applying the assumptions, description of terms above,

Table 1. Description of state variables of model (1)–(8).

Variables	Description
S	Low-risk susceptible individuals
J	High-risk susceptible individuals
E	Exposed individuals
A	Asymptomatic infectious individuals
B	Symptomatic infectious individuals
H	Hospitalized individuals
C	Death class
R	Recovered individuals

,

Table 2. Summary of the parameters.

Parameter	Meaning	Value	Reference
$\sigma$	Incubation rate	0.166	[15]
$\alpha$	Transmission rate between low-risk/high-risk susceptibles and asymptomatics	0.001	Assumed
$\beta$	Transmission rate between low-risk/high-risk susceptibles and symptomatics	0.001	Assumed
$\psi$	Fraction of low-risk/high-risk susceptibles who adhere to prevention guidelines	0.1	[22]
$\nu$	Fraction of low-risk/high-risk susceptibles who are on lock-down	0.7	Assumed
$ϵ$	Fraction of the exposed who become asymptomatic	0.2	[14]
m	Lock-down efficacy for low-risk susceptibles	0.5	[22]
g	Effectiveness of adherence for low-risk/high-risk susceptibles	0.5	[22]
p	Fraction of the exposed who are hospitalized	0.5	Assumed
q	Rate of testing	0.5	Assumed
$\tau$	Recovery rate for asymptomatics	$0.14286$	[14]
$\phi$	Recovery rate for hospitalized	$0.05$	[15]
d	Recovery rate for symptomatics	$0.1$	[15]
$\mathrm{\Lambda }$	Recruitment rate	10000	Assumed
$\varphi$	Vertical transmission recruitment rate	0.5	Assumed
$\mu$	Natural death rate of individuals	$0.9$	Assumed
$\delta$	COVID-19-caused death rate for symptomatics	0.015	[15]
u	COVID-19-caused death rate for hospitalized	0.015	[15]
k	Hospitalization rate for symptomatics	0.2	[15]
$\kappa$	Hospitalization rate for asymptomatics	0.19466	[15]
$\pi$	Reduction factor in COVID-19 transmission for concerned symptomatics	0.0242 (0–1)	[23,24]
a	Rate of safe burial of deaths	0.7	Assumed

and Figure 1, the differential equations which describe the dynamics of COVID-19 are formulated below.

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dS}{dt}}& =& \mathrm{\Lambda }-\alpha \pi (1-\nu m)(1-\psi g)SA-\beta (1-\nu m)(1-\psi g)SB-\mu S,\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dJ}{dt}}& =& \varphi -\alpha \pi (1-\nu m)(1-\psi g)JA-\beta (1-\nu m)(1-\psi g)JB-\mu J,\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dE}{dt}}& =& \alpha \pi (1-\nu m)(1-\psi g)SA+\beta (1-\nu m)(1-\psi g)SB+\hfill \\ & & \alpha \pi (1-\nu m)(1-\psi g)JA+\beta (1-\nu m)(1-\psi g)JB-(qp+\sigma (1-p)+\mu )E,\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dA}{dt}}& =& (1-p)ϵ\sigma E-(\kappa +\tau +\mu )A,\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dB}{dt}}& =& (1-p)(1-ϵ)\sigma E-(k+\delta +d+\mu )B,\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dH}{dt}}& =& qpE+\kappa A+kB-(u+\phi +\mu )H,\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dC}{dt}}& =& \delta B+uH-aC,\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dR}{dt}}& =& dB+\phi H+\tau A-\mu R.\hfill \end{array}$$

(8)

with initial conditions

$$\begin{array}{c}\hfill S\left(0\right)=S_0>0,J\left(0\right)=J_0>0,E\left(0\right)=E_0>0,A\left(0\right)=A_0>0,\\ \hfill B\left(0\right)=B_0>0,H\left(0\right)=H_0>0,C\left(0\right)=C_0,R\left(0\right)=R_0>0\end{array}$$

where

$$N\left(t\right)=S\left(t\right)+J\left(t\right)+E\left(t\right)+A\left(t\right)+B\left(t\right)+H\left(t\right)+C\left(t\right)+R\left(t\right)$$

2.1. Positivity of Solutions

Theorem 1.

Let $N\left(0\right)=N_0$ and let the initial data for (1)–(8) be $S\left(0\right)\ge 0,J\left(0\right)\ge 0,E\left(0\right)\ge 0,A\left(0\right)\ge 0,B\left(0\right)\ge 0,H\left(0\right)\ge 0,C\left(0\right)\ge 0$ and $R\left(0\right)\ge 0$. Then, the solutions $\left(S\right(t),J(t),E(t),A(t),B(t),H(t),C(t),R(t\left)\right)$ of the model will remain non-negative for all time $t>0$.

Proof. 

From (1), it can be seen that

${\displaystyle \frac{dS}{S}}=-(\alpha \pi (1-\nu m)(1-\psi g)A+\beta (1-\nu m)(1-\psi g)B+\mu )dt$

Further simplification gives

$$S\left(t\right)=S\left(0\right)exp(-{\int }_0^t(\alpha \pi (1-\nu m)(1-\psi g)A+\beta (1-\nu m)(1-\psi g)B+\mu )ds)>0\mathrm{for}t>0$$

Furthermore, from (2), it can be seen that

$${\displaystyle \frac{dJ}{J}}=-(\alpha \pi (1-\nu m)(1-\psi g)A+\beta (1-\nu m)(1-\psi g)B+\mu )dt$$

Further simplification gives

$J\left(t\right)=J\left(0\right)exp(-{\int }_0^t(\alpha \pi (1-\nu m)(1-\psi g)A+\beta (1-\nu m)(1-\psi g)B+\mu )>0$ for $t>0$

We can use similar arguments to show that $E\left(0\right)\ge >0,A\left(0\right)\ge 0,B\left(0\right)\ge 0,H\left(0\right)\ge 0,C\left(0\right)\ge 0,R\left(0\right)\ge 0$ □

Therefore, the model (1)–(8) is well posed and epidemiologically feasible in the region

$$D=\left((S,J,E,A,B,H,C,R)\in {\Re }_{+}^8:0\le S\le {\displaystyle \frac{\mathrm{\Lambda }}{\mu }},0\le J\le {\displaystyle \frac{\varphi }{\mu }},N\le {\displaystyle \frac{\mathrm{\Lambda }}{\mu }},S_0\le {\displaystyle \frac{\mathrm{\Lambda }}{\mu }},J_0\le {\displaystyle \frac{\varphi }{\mu }}\right)$$

Equations (1)–(5) are independent of H, C, and R, and therefore we decouple the compartments H, C, and R from the above equations which gives

Remark 1.

Models (1)–(8) are potentially analytically solvable by using the Lie algebra method through matrix exponentials. See in [25] and the references therein. Based on the Kolmogorov equation and the Wei–Norman method, the analytical solution for the proposed model can be obtained in terms of matrix exponentials.

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dS}{dt}}& =& \mathrm{\Lambda }-\alpha \pi (1-\nu m)(1-\psi g)SA-\beta (1-\nu m)(1-\psi g)SB-\mu S,\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dJ}{dt}}& =& \varphi -\alpha \pi (1-\nu m)(1-\psi g)JA-\beta (1-\nu m)(1-\psi g)JB-\mu J,\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dE}{dt}}& =& \alpha \pi (1-\nu m)(1-\psi g)SA+\beta (1-\nu m)(1-\psi g)SB+\hfill \\ \hfill & & \alpha \pi (1-\nu m)(1-\psi g)JA+\beta (1-\nu m)(1-\psi g)JB-(qp+\sigma (1-p)+\mu )E,\hfill \end{array}$$

(11)

$${\displaystyle \frac{dA}{dt}}=(1-p)ϵ\sigma E-(\kappa +\tau +\mu )A,$$

(12)

$${\displaystyle \frac{dB}{dt}}=(1-p)(1-ϵ)\sigma E-(k+\delta +d+\mu )B.$$

(13)

2.2. Basic Reproduction Number

First, we obtain the disease-free equilibrium point as $E_0=(S,J,E,A,B)=(S_0,J_0,0,0,0)=(\frac{\mathrm{\Lambda }}{\mu },\frac{\varphi }{\mu },0,0,0)$. Calculation of $E_0$ gives $\left({\displaystyle \frac{10,000}{0.9}},{\displaystyle \frac{0.5}{0.9}},0,0,0\right)$. This gives $(11111.11,0.555,0,0,0)$. We use the next-generation operator approach in [26] to obtain the basic reproduction number $R_0$. The approach is explored on the model (9)–(13) as follows:

Considering only the disease compartments, the rate of appearance of new infections and the transition rate are given, respectively, by

$$\left[\begin{array}{c}{\displaystyle \frac{dE}{dt}}\\ \\ {\displaystyle \frac{dA}{dt}}\\ \\ {\displaystyle \frac{dB}{dt}}\end{array}\right]=F-G=\left[\begin{array}{c}\alpha \pi (1-\nu m)(1-\psi g)SA+\beta (1-\nu m)(1-\psi g)SB+H_T\\ 0\\ 0\end{array}\right]-\left[\begin{array}{c}{k}_{1}E\\ -(1-p)ϵ\sigma E+k_{2}A\\ -(1-p)(1-ϵ)\sigma E+k_{3}B\end{array}\right]$$

where

$H_T=\alpha \pi (1-\nu m)(1-\psi g)JA+\beta (1-\nu m)(1-\psi g)JB$

$R_0$ is the maximum eigenvalue given as

$$\begin{array}{c}\hfill R_0^2={\displaystyle \frac{\alpha \pi (1-\nu m)(1-\psi g)(1-p)S_0ϵ\sigma }{k_1{k}_2}}+\end{array}$$

$$\begin{array}{c}\hfill {\displaystyle \frac{\alpha \pi (1-\nu m)(1-\psi g)(1-p)J_0ϵ\sigma }{k_1{k}_2}}+\end{array}$$

$$\begin{array}{c}\hfill {\displaystyle \frac{\beta (1-\nu m)(1-\psi g)(1-p)(1-ϵ)S_0\sigma }{k_1{k}_3}}+\end{array}$$

$$\begin{array}{c}\hfill {\displaystyle \frac{\beta (1-\nu m)(1-\psi g)(1-p)(1-ϵ)J_0\sigma }{k_1{k}_3}}\end{array}$$

This is the same as

$$\begin{array}{c}\hfill R_0^2={\displaystyle \frac{\alpha \pi (1-\nu m)(1-\psi g)(1-p)ϵ\sigma (S_0+J_0)}{k_1{k}_2}}+\end{array}$$

$$\begin{array}{c}\hfill {\displaystyle \frac{\beta (1-\nu m)(1-\psi g)(1-p)(1-ϵ)\sigma (S_0+J_0)}{k_1{k}_3}}\end{array}$$

Further simplification of $R_0$ gives

$$\begin{array}{c}\hfill R_0^2={\displaystyle \frac{\alpha \pi ϵ\sigma (1-\nu m)(1-\psi g)(1-p)(\mathrm{\Lambda }+\varphi )}{\mu (qp+\sigma (1-p)+\mu )(\kappa +\tau +\mu )}}+\end{array}$$

$$\begin{array}{c}\hfill {\displaystyle \frac{\beta \sigma (1-\nu m)(1-\psi g)(1-p)(1-ϵ)(\mathrm{\Lambda }+\varphi )}{\mu (qp+\sigma (1-p)+\mu )(\kappa +\delta +d+\mu )}}\end{array}$$

for

$$\begin{array}{ccc}\hfill k_1& =& qp+\sigma (1-p)+\mu \hfill \\ \hfill k_2& =& \kappa +\tau +\mu \hfill \\ \hfill k_3& =& k+\delta +d+\mu \hfill \end{array}$$

3. Global Stability Analysis

In reality, there are always small perturbations. It is of interest to know whether such deviations from steady state will lead to drastic changes or will be damped out in the advancement of time. Therefore, in what follows, we establish a global asymptotic stability of the disease-free equilibrium, that is, solution trajectories of sub-population, irrespective of where they start, converge to the disease-free equilibrium solution for $R_0<1$. This is stated and proved below.

3.1. Global Stability of Disease-Free Equilibrium Solution

Theorem 2.

The disease-free equilibrium solution is globally asymptotically stable if $R_0<1$ and unstable if $R_0>1$

Proof. 

Consider the following linear Lyapunov function

$$\begin{array}{c}\hfill L_1=E+{\displaystyle \frac{\alpha \pi (1-\nu m)(1-\psi g)(S_0+J_0)}{k_2}}A+{\displaystyle \frac{\beta (1-\nu m)(1-\psi g)(S_0+J_0)}{k_3}}B\end{array}$$

Its time derivative along solution path of COVID-19 model (9)–(13) is

$$\begin{array}{ccc}\hfill {\dot{L}}_1& =& \alpha \pi (1-\nu m)(1-\psi g)SA+\beta (1-\nu m)(1-\psi g)SB+\hfill \\ & & \alpha \pi (1-\nu m)(1-\psi g)JA+\beta (1-\nu m)(1-\psi g)JB-k_{1}E+\hfill \\ & & {\displaystyle \frac{\alpha \pi (1-\nu m)(1-\psi g)S_0}{k_2}}((1-p)ϵ\sigma E-k_{2}A)+\hfill \\ & & {\displaystyle \frac{\alpha \pi (1-\nu m)(1-\psi g)J_0}{k_2}}((1-p)ϵ\sigma E-k_{2}A)+\hfill \\ & & {\displaystyle \frac{\beta (1-\nu m)(1-\psi g)S_0}{k_3}}((1-p)(1-ϵ)\sigma E-k_{3}B)+\hfill \\ & & {\displaystyle \frac{\beta (1-\nu m)(1-\psi g)J_0}{k_3}}((1-p)(1-ϵ)\sigma E-k_{3}B)\hfill \\ \hfill & =& S\left(\alpha \pi \right(1-\nu m\left)\right(1-\psi g)A+\beta (1-\nu m\left)\right(1-\psi g)B-\hfill \\ & & S_0(\alpha \pi (1-\nu m)(1-\psi g)A+\beta (1-\nu m)(1-\psi g)B+\hfill \\ & & J\left(\alpha \pi \right(1-\nu m\left)\right(1-\psi g)A+\beta (1-\nu m\left)\right(1-\psi g\left)B\right)-\hfill \\ & & J_0(\alpha \pi (1-\nu m)(1-\psi g)A+\beta (1-\nu m)(1-\psi g)B)\hfill \\ & & +\left({\displaystyle \frac{D_T}{k_1{k}_2}}+{\displaystyle \frac{F_T}{k_1{k}_2}}+{\displaystyle \frac{W_T}{k_1{k}_3}}+{\displaystyle \frac{N_T}{k_1{k}_3}}-1\right)k_{1}E\hfill \\ \hfill & =& (S-S_0)(\alpha \pi (1-\nu m)(1-\psi g)A+\beta (1-\nu m)(1-\psi g)B)+\hfill \\ & & (J-J_0)(\alpha \pi (1-\nu m)(1-\psi g)A+\beta (1-\nu m)(1-\psi g)B)+(R_0^2-1)k_{1}E\hfill \end{array}$$

where

$$\begin{array}{ccc}\hfill D_T& =& S_0ϵ\sigma \alpha \pi (1-\nu m)(1-\psi g)(1-p)\hfill \\ \hfill F_T& =& J_0\alpha \pi (1-\nu m)(1-\psi g)(1-p)ϵ\sigma \hfill \\ \hfill W_T& =& S_0\beta \sigma (1-\nu m)(1-\psi g)(1-p)(1-ϵ)\hfill \\ \hfill N_T& =& J_0\beta \sigma (1-\nu m)(1-\psi g)(1-p)(1-ϵ)\hfill \end{array}$$

Because $S\le S_0$ and $J\le J_0$, we have

$${\dot{L}}_1=k_{1}E(R_0^2-1)\le 0$$

whenever $R_0\le 1$

Therefore, ${\dot{L}}_1\le 0$ for $R_0\le 1$ and ${\dot{L}}_1=0$ if and only if $E=0,A=0,B=0$ or $S=S_0$, $J=J_0$ and $R_0=1$. Consequently, the largest invariant set in $(S,J,E,A,B,H,C,R)\in D:L_1^{\prime }=0$ is the singleton $E_0$ and by LaSalle’s Invariance Principle [27], $E_0$ is globally asymptotically stable. The epidemiological implication of the above result is that COVID-19 can be eradicated from the population when $R_0\le 1$ irrespective of the initial sizes of the sub-population of the model. □

3.2. Global Stability of Endemic Equilibrium Point

Theorem 3.

The model (9)–(13) has a unique endemic equilibrium solution whenever $R_0>1$, and no endemic equilibrium solution otherwise.

Proof. 

The following equations are satisfied by the endemic equilibrium solution $E_1=(S^{*},J^{*},E^{*},A^{*},B^{*})$:

$$\mathrm{\Lambda }-\alpha \pi (1-\nu m)(1-\psi g)S^{*}{A}^{*}-\beta (1-\nu m)(1-\psi g)S^{*}{B}^{*}-\mu S^{*}=0,$$

(14)

$$\varphi -\alpha \pi (1-\nu m)(1-\psi g)J^{*}{A}^{*}-\beta (1-\nu m)(1-\psi g)J^{*}{B}^{*}-\mu J^{*}=0,$$

(15)

$$\alpha \pi (1-\nu m)(1-\psi g)S^{*}{A}^{*}+\beta (1-\nu m)(1-\psi g)S^{*}{B}^{*}$$

$$+\alpha \pi (1-\nu m)(1-\psi g)J^{*}{A}^{*}+\beta (1-\nu m)(1-\psi g)J^{*}{B}^{*}-k_1{E}^{*}=0$$

(16)

$$(1-p)ϵ\sigma E^{*}-k_2{A}^{*}=0,$$

(17)

$$(1-p)(1-ϵ)\sigma E^{*}-k_3{B}^{*}=0.$$

(18)

Making $S^{*}$, $J^{*}$, $A^{*}$, and $B^{*}$ the subject of formula in (14)–(18) gives

$$S^{*}={\displaystyle \frac{\mathrm{\Lambda }}{\alpha \pi (1-\nu m)(1-\psi g)A^{*}+\beta (1-\nu m)(1-\psi g)B^{*}+\mu }}$$

(19)

$$J^{*}={\displaystyle \frac{\varphi }{\alpha \pi (1-\nu m)(1-\psi g)A^{*}+\beta (1-\nu m)(1-\psi g)B^{*}+\mu }}$$

(20)

$$A^{*}={\displaystyle \frac{ϵ\sigma (1-p)E^{*}}{k_2}}$$

(21)

$$B^{*}={\displaystyle \frac{(1-p)(1-ϵ)\sigma E^{*}}{k_3}}$$

(22)

Adding (14)–(16), we get

$$\mathrm{\Lambda }+\varphi -\mu S^{*}-\mu J^{*}-k_1{E}^{*}=0$$

(23)

Substituting (19)–(22) in (23) yields

$$-{\displaystyle \frac{k_1\mu R_0^2{\left(E^{*}\right)}^2}{\mathrm{\Lambda }+\varphi }}+\mu (R_0^2-1)E^{*}=0$$

(24)

From (24), we can obtain $E^{*}$ to be

$$E^{*}={\displaystyle \frac{\mathrm{\Lambda }+\varphi }{k_1}}\left(1-{\displaystyle \frac{1}{R_0^2}}\right)$$

Substituting $E^{*}$ in (19)–(22) gives

$$S^{*}={\displaystyle \frac{\mathrm{\Lambda }}{\mu R_0^2}}$$

(25)

$$J^{*}={\displaystyle \frac{\varphi }{\mu R_0^2}}$$

(26)

$$A^{*}={\displaystyle \frac{(1-p)(\mathrm{\Lambda }+\varphi )ϵ\sigma }{k_1{k}_2}}\left(1-{\displaystyle \frac{1}{R_0^2}}\right)$$

(27)

$$B^{*}={\displaystyle \frac{(1-p)(1-ϵ)(\mathrm{\Lambda }+\varphi )\sigma }{k_1{k}_3}}\left(1-{\displaystyle \frac{1}{R_0^2}}\right)$$

(28)

Calculation of $E^{*}$ yields

$$\begin{array}{ccc}\hfill S^{*}& =& 41534.09533\hfill \\ \hfill J^{*}& =& 0.06922349223\hfill \\ \hfill A^{*}& =& 2857.065333\hfill \\ \hfill B^{*}& =& 11640.08393\hfill \end{array}$$

The above results exist if only $R_0>1$.

Further, we establish the global asymptotic stability of the endemic equilibrium solution by proving that solution trajectories converge to the endemic equilibrium point for $R_0>1$. We shall carry this out by constructing a suitable Lyapunov candidate function of Goh–Volterra type (see [28,29] and the references cited therein). The result below establishes the global stability of the endemic equilibrium solution $E_1$. □

Theorem 4.

If $R_0>1$, then the endemic equilibrium $E_1$ of the model (9)–(13) is globally asymptotically stable in the interior of the region D.

Proof. 

Given the following equations which are satisfied by the endemic equilibrium point $E_1$:

$$\mathrm{\Lambda }=\alpha \pi (1-\nu m)(1-\psi g)S^{*}{A}^{*}+\beta (1-\nu m)(1-\psi g)S^{*}{B}^{*}+\mu S^{*},$$

(29)

$$\varphi =\alpha \pi (1-\nu m)(1-\psi g)J^{*}{A}^{*}+\beta (1-\nu m)(1-\psi g)J^{*}{B}^{*}+\mu J^{*}$$

(30)

$$\alpha \pi (1-\nu m)(1-\psi g)S^{*}{A}^{*}+\beta (1-\nu m)(1-\psi g)S^{*}{B}^{*}$$

$$+\alpha \pi (1-\nu m)(1-\psi g)J^{*}{A}^{*}+\beta (1-\nu m)(1-\psi g)J^{*}{B}^{*}=k_1{E}^{*}$$

(31)

$$(1-p)ϵ\sigma E^{*}=k_2{A}^{*},$$

(32)

$$(1-p)(1-ϵ)\sigma E^{*}=k_3{B}^{*}.$$

(33)

Consider the following Lotka–Volterra Lyapunov function

$$\begin{array}{c}\hfill L_2=\left(S-S^{*}-S^{*}ln{\displaystyle \frac{S^{*}}{S^{*}}}\right)+\left(J-J^{*}-J^{*}ln{\displaystyle \frac{J^{*}}{J^{*}}}\right)+\left(E-E^{*}-E^{*}ln{\displaystyle \frac{E^{*}}{E^{*}}}\right)+\\ \hfill n_1\left(A-A^{*}-A^{*}ln{\displaystyle \frac{A}{A^{*}}}\right)+n_2\left(B-B^{*}-B^{*}ln{\displaystyle \frac{B^{*}}{B^{*}}}\right)\end{array}$$

where

$$n_1={\displaystyle \frac{\alpha \pi (1-\nu m)(1-\psi g)S^{*}}{k_2}}+{\displaystyle \frac{\alpha \pi (1-\nu m)(1-\psi g)J^{*}}{k_2}}$$

$$n_2={\displaystyle \frac{\beta (1-\nu m)(1-\psi g)S^{*}}{k_3}}+{\displaystyle \frac{\beta (1-\nu m)(1-\psi g)J^{*}}{k_3}}$$

The Lyapunov time derivative along the solutions of the model gives

$${\dot{L}}_2=\left(1-{\displaystyle \frac{S^{*}}{S}}\right)S^{\prime }+\left(1-{\displaystyle \frac{J^{*}}{J}}\right)J^{\prime }+\left(1-{\displaystyle \frac{E^{*}}{E}}\right)E^{\prime }+n_1\left(1-{\displaystyle \frac{A^{*}}{A}}\right)A^{\prime }+n_2\left(1-{\displaystyle \frac{B^{*}}{B}}\right)B^{\prime }$$

$$\begin{gathered} {\dot{L}}_2=\left(1-{\displaystyle \frac{S^{*}}{S}}\right)((\alpha \pi (1-\nu m)(1-\psi g)A^{*}+\beta (1-\nu m)(1-\psi g)B^{*})S^{*}+\mu S^{*}-(\alpha \pi (1- \\ \nu m)(1-\psi g)A+\beta (1-\nu m)(1-\psi g)B)S-\mu S+\left(1-{\displaystyle \frac{J^{*}}{J}}\right)((\alpha \pi (1-\nu m)(1-\psi g)A^{*}+\beta (1- \\ \nu m)(1-\psi g)B^{*})J^{*}+\mu J^{*}-(\alpha \pi (1-\nu m)(1-\psi g)A+\beta (1-\nu m)(1-\psi g)B)J-\mu J+ \\ \left(1-{\displaystyle \frac{E^{*}}{E}}\right)((\alpha \pi (1-\nu m)(1-\psi g)SA+\beta (1-\nu m)(1-\psi g)SB+\alpha \pi (1-\nu m)(1-\psi g)JA+ \\ \beta (1-\nu m)(1-\psi g)JB-k_{1}E)+n_1\left(1-{\displaystyle \frac{A^{*}}{A}}\right)((1-p)ϵ\sigma E-k_{2}A)+n_2\left(1-{\displaystyle \frac{B^{*}}{B}}\right)((1-p)(1- \\ ϵ)\sigma E-k_{3}B) \end{gathered}$$

Further simplification gives

$$\begin{gathered} {\dot{L}}_2=(\alpha \pi (1-\nu m)(1-\psi g)A^{*}+\beta (1-\nu m)(1-\psi g)B^{*})S^{*} \\ +(\alpha \pi (1-\nu m)(1-\psi g)A^{*}+\beta (1-\nu m)(1-\psi g)B^{*})J^{*}+k_1{E}^{*}+n_1{k}_2{A}^{*}+n_2{k}_3{B}^{*}+p_3{k}_4{U}^{*} \\ -(\alpha \pi (1-\nu m)(1-\psi g)A^{*}+\beta (1-\nu m)(1-\psi g)B^{*}{\displaystyle \frac{{\left(S^{*}\right)}^2}{S}}-(\alpha \pi (1-\nu m)(1-\psi g)A^{*}+\beta (1- \\ \nu m)(1-\psi g)B^{*}{\displaystyle \frac{{\left(J^{*}\right)}^2}{J}}-(\alpha \pi (1-\nu m)(1-\psi g)A+\beta (1-\nu m)(1-\psi g)B{\displaystyle \frac{S{E}^{*}}{E}}-(\alpha \pi (1-\nu m) \\ (1-\psi g)A+\beta (1-\nu m)(1-\psi g)B{\displaystyle \frac{J{E}^{*}}{E}}-{\displaystyle \frac{n_1(1-p)ϵ\sigma E{A}^{*}}{A}}-{\displaystyle \frac{n_2(1-p)(1-ϵ)\sigma E{B}^{*}}{B}}+2\mu S^{*} \\ +2\mu J^{*}-{\displaystyle \frac{\mu {\left(S^{*}\right)}^2}{S}}-{\displaystyle \frac{\mu {\left(J^{*}\right)}^2}{J}}-\mu (S+J)+(\alpha \pi (1-\nu m)(1-\psi g)A+\beta (1-\nu m)(1-\psi g)B \\ +\mu )S^{*}+(\alpha \pi (1-\nu m)(1-\psi g)A+\beta (1-\nu m)(1-\psi g)B+\mu )J^{*} \end{gathered}$$

Substituting $n_1$ and $n_2$ by their values and employing (29)–(33) gives

$$n_1(1-p)ϵ\sigma ={\displaystyle \frac{\alpha \pi (1-\nu m)(1-\psi g)S^{*}{A}^{*}+\alpha \pi (1-\nu m)(1-\psi g)J^{*}{A}^{*}}{E^{*}}}$$

(34)

$$n_2(1-p)(1-ϵ)\sigma ={\displaystyle \frac{\beta (1-\nu m)(1-\psi g)S^{*}{B}^{*}+\beta (1-\nu m)(1-\psi g)J^{*}{B}^{*}}{E^{*}}}$$

(35)

Using (29)–(35), we have

$$\begin{gathered} {\dot{L}}_2=\mu S^{*}\left(2-{\displaystyle \frac{S^{*}}{S}}+\mu J^{*}\left(2-{\displaystyle \frac{J^{*}}{J}}-{\displaystyle \frac{J}{J^{*}}}\right)\right)+3(\alpha \pi (1-\nu m)(1-\psi g)A^{*}+\beta (1-\nu m) \\ (1-\psi g)B^{*})S^{*}+3(\alpha \pi (1-\nu m)(1-\psi g)A^{*}+\beta (1-\nu m)(1-\psi g)B^{*})J^{*}-(\alpha \pi (1-\nu m) \\ (1-\psi g)A^{*}+\beta (1-\nu m)(1-\psi g)B^{*}) \\ {\displaystyle \frac{{\left(S^{*}\right)}^2}{S}}-(\alpha \pi (1-\nu m)(1-\psi g)A^{*}+\beta (1-\nu m)(1-\psi g)B^{*}){\displaystyle \frac{{\left(J^{*}\right)}^2}{J}}-(\alpha \pi (1-\nu m)(1-\psi g) \\ A+\beta (1-\nu m)(1-\psi g)B){\displaystyle \frac{S{E}^{*}}{E}}-(\alpha \pi (1-\nu m)(1-\psi g)A+\beta (1-\nu m)(1-\psi g)B) \\ {\displaystyle \frac{J{E}^{*}}{E}}-{\displaystyle \frac{\alpha \pi (1-\nu m)(1-\psi g)S^{*}{E}^{*}{\left(A^{*}\right)}^2}{E^{*}A}} \\ -{\displaystyle \frac{\alpha \pi (1-\nu m)(1-\psi g)J^{*}{E}^{*}{\left(A^{*}\right)}^2}{E^{*}A}}- \\ {\displaystyle \frac{\beta (1-\nu m)(1-\psi g)S^{*}E{\left(B^{*}\right)}^2}{E^{*}B}}-{\displaystyle \frac{\beta (1-\nu m)(1-\psi g)J^{*}E{\left(B^{*}\right)}^2}{E^{*}B}} \end{gathered}$$

$$\begin{gathered} {\dot{L}}_2=\mu S^{*}\left(2-{\displaystyle \frac{S^{*}}{S}}-{\displaystyle \frac{S}{S^{*}}}\right)+\mu J^{*}\left(2-{\displaystyle \frac{J^{*}}{J}}-{\displaystyle \frac{J}{J^{*}}}\right)+ \\ \alpha \pi (1-\nu m)(1-\psi g)S^{*}{A}^{*}\left(3-{\displaystyle \frac{S^{*}}{S}}-{\displaystyle \frac{AS{E}^{*}}{A^{*}E{S}^{*}}}-{\displaystyle \frac{A^{*}{E}^{*}}{E^{*}A}}\right)+ \\ \alpha \pi (1-\nu m)(1-\psi g)J^{*}{A}^{*}\left(3-{\displaystyle \frac{J^{*}}{J}}-{\displaystyle \frac{AJ{E}^{*}}{A^{*}E{J}^{*}}}-{\displaystyle \frac{A^{*}{E}^{*}}{E^{*}A}}\right)+ \\ \beta (1-\nu m)(1-\psi g)S^{*}{B}^{*}\left(3-{\displaystyle \frac{S^{*}}{S}}-{\displaystyle \frac{BS{E}^{*}}{B^{*}E{S}^{*}}}-\right.{\displaystyle \frac{E{B}^{*}}{E^{*}B}}+\beta (1-\nu m)(1-\psi g)J^{*}{B}^{*} \\ \left(3-{\displaystyle \frac{J^{*}}{J}}-{\displaystyle \frac{BJ{E}^{*}}{B^{*}E{J}^{*}}}-{\displaystyle \frac{E{B}^{*}}{E^{*}B}}\right) \end{gathered}$$

By arithmetic-geometric means inequality, i.e., $n-(a_1+a_2+\cdots +a_n)\le 0$, where $a_1.a_2\dots a_n=1$ and $a_1,a_2,\dots ,a_n>0$, it follows that ${\dot{L}}_2\le 0$ with $L_2=0$ if and only if $S=S^{*}$, $J=J^{*}$, $E=E^{*}$, $A=A^{*}$, $B=B^{*}$.

Therefore, the largest compact invariant subset of the set where ${\dot{L}}_2=0$ is

$$(S,J,E,A,B)=(S^{*},J^{*},E^{*},A^{*},B^{*})$$

Therefore, by stability theorem of Lyapunov and LaSalle’s Invariance Principle, it follows that every solution in D approaches $E_1$ for $R_0>1$ as $t\to \infty$.

The implication of the above result is that COVID-19 will establish itself in the population for $R_0>1$. □

4. Sensitivity Analysis

We perform a sensitivity analysis of COVID-19 model in order to determine the parameters that are most influential in the transmission dynamics and spread of the disease. Sensitivity analysis is helpful for experimental design, data assimilation, and reduction of complex nonlinear models. A very high sensitivity indicates that more care should be taken in the estimation of the associated parameters. Following the approach in [30], we define sensitivity index of a variable to a parameter as the ratio of relative change in the variable to the relative change in the parameter.

When the variable is a differentiable function of the parameter, the sensitivity index may be alternatively defined using partial derivatives.

Definition 1.

The normalized forward sensitivity index of a variable, u(p), that depends differentiably on a parameter, p, is defined as

$$N_p^u=\frac{\partial u}{\partial p}\times \frac{p}{u}$$

for $u\ne 0$

Consequently, we derive analytical expression for the sensitivity index of $R_0$ as

$$N_{p_i}^{R_0}=\frac{\partial R_0}{\partial p_i}\times \frac{p_i}{R_0}$$

where $p_i,i\in \mathbb{N}$ denotes each parameter involved in $R_0$

$R_0$ is defined as

$$\begin{array}{c}\hfill R_0^2={\displaystyle \frac{\alpha \pi ϵ\sigma (1-\nu m)(1-\psi g)(1-p)(\mathrm{\Lambda }+\varphi )}{\mu (qp+\sigma (1-p)+\mu )(\kappa +\tau +\mu )}}+\end{array}$$

$$\begin{array}{c}\hfill {\displaystyle \frac{\beta \sigma (1-\nu m)(1-\psi g)(1-p)(1-ϵ)(\mathrm{\Lambda }+\varphi )}{\mu (qp+\sigma (1-p)+\mu )(\kappa +\delta +d+\mu )}}\end{array}$$

We have

Table 3. Numerical values of sensitivity indices of $R_0$ with respect to parameter involved.

Parameter Symbol	Sensitivity Index
$\mathrm{\Lambda }$	+0.4999
$\beta$	+0.4963
$\sigma$	+0.4663
$\alpha$	+0.0036
$\varphi$	+0.00000083
$\pi$	+0.0036
$\mu$	−1.2352
p	−0.5689
m	−0.3321
$\nu$	−0.3308
$ϵ$	−0.1203
q	−0.1013
k	−0.0816
g	−0.0555
d	−0.0408
$\psi$	−0.0555
$\kappa$	−0.0005
$\delta$	−0.0061
$\tau$	−0.0004

which summarizes the sensitivity indices on $R_0$ with respect to parameters, i.e.,

$$\begin{array}{c}\hfill N_{\mathrm{\Lambda }}^{R_0},N_{\beta }^{R_0},N_{\sigma }^{R_0},N_{\alpha }^{R_0},N_{\varphi }^{R_0},N_{\pi }^{R_0},N_{\mu }^{R_0},N_p^{R_0},N_m^{R_0},\\ \hfill N_{\nu }^{R_0},N_{ϵ}^{R_0},N_q^{R_0},N_k^{R_0},N_g^{R_0},N_d^{R_0},N_{\psi }^{R_0},N_{\kappa }^{R_0},N_{\delta }^{R_0},N_{\tau }^{R_0}\end{array}$$

In order to know the most influential parameters concerned with the basic reproduction of the disease, sensitivity analysis is performed. We categorize the most sensitive parameters in descending order to the least sensitive ones in Table 3, after ranking the sensitivity functions. The most sensitive parameters that have positive indices on $R_0$ are the recruitment rate $\mathrm{\Lambda }$, followed by the transmission rate between low-risk/high-risk susceptible individuals and symptomatic infectious individuals $\beta$, followed by the incubation rate $\sigma$, followed by the transmission rate between low-risk/high-risk susceptible individuals and asymptomatic infectious individuals $\alpha$, followed by vertically transmitted recruitment rate $\varphi$. From the sensitivity analysis results, we can confidently say that the intervention strategies must focus on the reduction in the above parameters in order to have a disease-free population. On the other hand, the most sensitive parameters in a negative sense is the natural death rate $\mu$, followed by the fraction of the exposed individuals who are tested for COVID-19 p and become hospitalized, followed by the lockdown efficacy for low and high risk susceptible individuals m, followed by the fraction of low-risk and high-risk susceptible individuals who are on lockdown $\nu$, followed by the fraction of the exposed individuals who become asymptomatic $ϵ$, followed by rate of testing q, followed by hospitalization rate of symptomatic k, followed by effectiveness of adherence for low-risk and high-risk susceptible individuals g, followed by the recovery rate of symptomatic individuals d, followed by the fraction of low-risk and high-risk susceptible individuals who adhere to prevention guidelines $\psi$, followed by hospitalization rate of asymptomatic infectious individuals $\kappa$, followed by COVID-19-caused death rate $\delta$, followed by recovery rate of asymptomatic infectious individuals $\tau$. It is clear to see that a fraction p, of the exposed individuals who are tested for COVID-19 at a rate q, and become hospitalized, has significant contribution to the disease dynamics. Therefore, high level of contact tracing through testing and hospitalization of the exposed individuals who are tested for COVID-19, are needed at the beginning of the epidemic. The sensitivity indices of the model (9)–(12) are graphically shown in Figure 2.

5. Numerical Simulations

Numerical simulations for COVID-19 model are presented in this section with the parameters mostly taken from the literature. See Table 1 and Table 2 for the description of variables and parameters with their values. Figure 3 represents the solution trajectories of the model for initial conditions $\left(S\right(0)=50,J(0)=40,E(0)=10,A(0)=150,B\left(0\right)=200,$ $H\left(0\right)=100,$ $C\left(0\right)=10)$, $R\left(0\right)=150)$ when $\mathrm{\Lambda }=10,000$, $\varphi =0.5$ and $\mu =0.9$ so that $R_0=0.517$. From Figure 3, it can be seen that the population of exposed individuals E, asymptomatic infectious individuals A, symptomatic infectious individuals B, hospitalized individuals H, and death class C converge to disease-free equilibrium solution and never turn to epidemic for $R_0<1$ and this gives an increase in the population of low-risk susceptible individuals S and high- risk susceptible individuals. This implies that COVID-19 is cleared from the population for $R_0<1$ as shown in Theorem 2. It can be viewed that the graph of high risk susceptible individuals J do not rise as much as the graph of low risk susceptible individuals S. This may be due to the fact that the high risk susceptible individuals may prove a bit more difficult to control than the low risk susceptible individuals due to their constant and inevitable exposure to the virus. Figure 4 presents the solution trajectories when $\mathrm{\Lambda }$ = 300,000, $\varphi =0.5$ and $\mu =0.9$ for initial conditions $\left(S\right(0)=1000,000,J(0)=400,000,E(0)=500,$$A\left(0\right)=200,$$B\left(0\right)=250,H\left(0\right)=150,C\left(0\right)=20$, $R\left(0\right)=200)$ so that $R_0=2.832$. From Figure 4, it can be seen that irrespective of the initial size of the infectious population, the infectious individuals remain in the population for $R_0>1$. Therefore, the limit of the solution trajectories is the endemic equilibrium solution. This implies that solution trajectories converge to the endemic equilibrium point. Thus, as established in Theorem 4, the disease establishes itself in the population.

6. Conclusions

We have developed an epidemic model for the transmission dynamics of COVID-19. Previous research in the literature has concentrated on mathematical modeling of COVID-19, taking into account, the general public (susceptible population) without considering individuals with vertically transmitted underlying diseases (high-risk susceptible individuals), on the disease dynamics. Therefore, we formulate a mathematical model that incorporates a very important factor of individuals with higher risk of contacting the disease rather than the general public. These include HIV patients, HCV patients, asthma patient etc. The epidemiological threshold quantity, $R_0$, is derived and we establish that our model has a globally asymptotic stable disease-free equilibrium solution whenever $R_0<1$. For the case where $R_0>1$, we have proved that there exists a unique endemic equilibrium solution, which is proved analytically to be globally asymptotically stable. We illustrate that COVID-19 disappears from the population when $R_0<1$ but persists in the population for $R_0>1$. See Figure 3 and Figure 4. Numerical simulations reveal two equilibria: disease-free equilibrium solution and endemic equilibrium solution, which are shown to be globally asymptotically stable for $R_0<1$ and $R_0>1$. We perform the sensitivity analysis of the basic reproduction number to each of the parameters to explore the significance of each parameter to $R_0$. Sensitivity analysis results reveal that increase in the fraction of exposed individuals who are tested for COVID-19, through contact tracing, largely has negative impact on $R_0$. Therefore, in order to reduce the continuous increase in the spread of the disease in the population, increase in the effective contact tracing of the exposed individuals through testing must be activated. Then, the problem of limited hospital spaces should be addressed in order to absorb higher number of COVID-19 cases. It is also pertinent that the front-line health care workers who have underlying diseases, some of which are funeral directors and mortuary personnel, who handle the bodies of those confirmed to have died from COVID-19, should avoid unnecessary manipulation of the body that may expel air from the lungs. This study can help epidemiologists and public health officials to focus more on the parameters in formulating a disease control policy. This work can be extended by placing control on the most sensitive parameter in the model so as to know the duration of time to apply control for the disease to die out. This can be achieved through optimal control analysis. It can also be extended by considering the role of high risk susceptible population on the spread of COVID-19 in a meta-population model.