A Mathematical Model for Covid-19 Disease Transmission Dynamics with Impact of Saturated Treatment: Modeling, Analysis and Simulation

Idisi I. Oke¹, Yakub T. Oyebo¹, Olamide F. Fakoya¹, Vigbe S. Benson², Yusuf T. Tunde^3*

¹Department of Mathematics, Lagos State University, Lagos, Nigeria.
²Department of Technical Education, Alvan Ikoku Federal College of Education, Owerri, Imo State, Nigeria.
³Department of Mathematical Sciences, Federal University of Technology Akure, Akure, Nigeria.

DOI: 10.4236/oalib.1107332                

Abstract

The emergence of novel coronavirus (Covid-19) in 2019 sprung a sudden outbreak across the globe, presenting clinical and public health management challenges which led to global cancellation of conferences, travel restrictions, social distancing and closure of institutions. Thus, in considering the grave implications of the continuous spread of coronavirus disease, a SEIHRD epidemic model was formulated to gain insight into disease transmission dynamics with impacts of proposing control measures. The model captures the impact of undetected infectious individuals and detected hospitalized individuals with saturated treatment on the spread, death and recovery of Covid-19 patients in Nigeria. The model epidemic threshold and equilibria are determined, and their stabilities are analyzed. The model is validated by fitting it to data from January 28 to December 5, 2020. Results obtained suggest that decreasing the transmission rate for infective alone is not sufficient to eradicate the disease because of the presence of backward bifurcation, and we recommend that Nigerians must also adhere strictly to Covid-19 protocols in mitigating the spread and demise of coronavirus.

Keywords

Coronavirus, Basic Reproduction Number, Equilibrium Solutions, Stability Analysis, Backward Bifurcation

1. Introduction

The emergence of strain severe acute respiratory syndrome Coronavirus-2 (SARS-COV-2) commonly known as Covid-19 on December 31, 2019, sprung a sudden outbreak that span globally across (Europe, North America, Asia, South America, Africa and Oceania) affecting 215 countries/territories which rapidly becomes a devastating pandemic. As of June 2020, over 6,861,727 people were affected by coronavirus and over 3,361,466 deaths globally due to the virus [1]. The demise of Covid-19 disease is still on increase in bridging global health and social economy security, on the account of its devastating nature world health organization (WHO) on February 11, 2020, declared coronavirus to be Public Health Emergency of International Concern [2].

The outbreak of Covid-19 was reported in Wuhan, Hubei province, China on December 31, 2019, followed by Europe, America, Asia, Africa and across the globe. Coronavirus spreads rapidly across 47 African countries with Nigeria being the third largest affected African country having recorded over 84,811 confirmed cases and 1264 deaths. Nigeria first documented Covid-19 case was recorded on February 29^th 2020, after a 44-year-old Italian citizen on January 27^th 2020 arrived at the Murtala Muhammed International Airport, from Milan, Italy traveling to his company site at Ogun State. Starting with the official index from February, Nigerian began to experience increase in number of new cases, which led Nigeria government to respond to the outbreak by adhering strictly to Nigeria Centre for Disease Control (NCDC) and World Health Organization (WHO) preventive guidelines via enforcing total lock-down on human and economy activities in states with high index cases, while imposing travel restrictions, social distancing, closure of institutions and others unprecedented preventive measures in order to mitigate the spread of disease.

Coronaviruses are large family of zoonotic virus that causes disease in both humans and animals, In humans, there are three existing families known as coronaviruses which includes severe acute Respiratory Syndrome SARS-CoV discovered in (2002), Middle Eastern Respiratory Syndrome MERS-CoV (2012) and Covid-19 caused by SARS-CoV-2 (2019). In animals like cows and pigs, they cause diarrhea, while in mice they cause hepatitis and encephalomyelitis. The history of these viruses have their members generally associated with vector-transmitted diseases in humans and causes illness ranging from the common cold to severe respiratory disease which includes symptoms such as dry cough, fever, shortness of breath, and breathing difficulties. However, in more severe cases infection can cause pneumonia, kidney failure, SARS and death [3]. Coronavirus infection is transmitted from one person to another through droplets produced from the respiratory system of infected persons, often during coughing or sneezing and through inhalation of respiratory droplets from symptomatic and infectious humans. There is also limited evidence that the virus can be exhaled through normal breathing. The incubation period of Covid-19 ranges from 2 - 14 days, and most infections (over 80%) show mild or no clinical symptoms of the disease [4]. The ongoing outbreak presents many clinical and public health management challenges due to limited understanding of viral pathogen, risk factors for infection, period of infectivity, modes and extent of virus inter-human transmission, effective preventive measures, limited public health response and containment interventions coupled with no antiviral treatment nor vaccine available before the first authorization of the Pfizer-Bio NTech vaccine on December 11, 2020, which followed by the Moderna vaccine authorized by Fedral Drug Administration (FDA) on December 18, 2020, and these vaccines are yet to be duly circulated in Europe much more reaching Africa precisely Nigeria.

Managing this serious epidemic requires the appropriate deployment of limited human resources across all cadres of health care and public health as recommended in [5] Due to the un-adhesiveness of preventive measures and unpreparedness of some countries, the novel coronavirus has invaded some countries more than others and this has imperatively led to a high increase of infected cases. However, a variety of mathematical models [4] [6] [7] [8] [9] [10] have been developed to further gain insight in understanding the impact and spread of Covid-19 pandemic, some of these models such as Khan and Atangana [8], formulated a mathematical model to study the dynamics of novel coronavirus with fractional derivative using available data of infectious cases from Jan 21, 2020 to Jan 28, 2020 in Wuhan China, to parameterized their model. Ivora, B., et al. [10] developed a new epidemic model to study the dynamic of Coronavirus disease which captures the existence of undetected infectious cases and the different sanitary and infectiousness conditions of hospitalized people. Qianying et al. [9] formulated a conceptual model for the coronavirus disease (Covid-19) outbreak in Wuhan China, with individual reaction and governmental action to proffer insight to understanding the trend of coronavirus disease outbreak. While Iboi et al. [7], formulated a mathematical model to propel insight to Covid-19 pandemic in Nigeria and they concluded that Nigeria, should ramp up widespread diagnostic Covid-19 testing and contact tracing to have a realistic measure of the burden of the pandemic with emphasis on personal hygiene, hand washing (using soaps and approved hand sanitizers), physical-distancing, wearing face masks in public and make personal protection equipment widely available to frontline healthcare workers.

However, each country or city has its own shortfall on precautionary measures and maximal capacity for the treatment of coronavirus disease. This study is aimed to rigorously evaluate the impact of undetected infectious patients and limited supplies of medical facilities for treatment of hospitalized or isolated Covid-19 patients in Nigeria. Hence, this paper is organized as follows. The proposed model is formulated in Section (2) and Section (3) for model properties and analysis. Section (4) describes the steady-state analysis, Section (5) describes the model fitting and sensitivity analysis, Section (6) describes the numerical analysis and the conclusion is presented in Section (7).

2. Model Formulation

To effectively evaluate Covid-19 disease dynamics with saturated treatment, a deterministic mathematical model (SEIHRD) comprising six distinct compartments of susceptible $S\left(t\right)$, exposed $E\left(t\right)$, infected $I\left(t\right)$, hospitalized $H\left(t\right)$, recovered $R\left(t\right)$ and deceased $D\left(t\right)$ is formulated with inclusion of saturated treatment function as applicable to similar pandemic diseases see [11] [12]. Susceptible individuals are recruited into the population at rate A, and are exposed to Covid-19 virus at (rate of $\lambda$ ). Susceptible became infected with the pathogen with no clinical symptoms while undergoing disease incubation circle, an (average of 5-days), after which the exposed individual becomes infectious and can transmit coronavirus. The infected person either maintained self isolation or notified the medical response authority (team) to be hospitalized with access to Covid-19 treatment facilities. The self isolated individual enters into the compartment $I\left(t\right)$ at a transformation (rate a), who in some rear cases recovered from the disease at (rate ${\omega }_1$ ) or die from Covid-19 at (rate ${\delta }_1$ ). In this case, these individuals are assume to be highly infective and are undetected by the health authority. Hospitalized Covid-19 patients enters into compartment $H\left(t\right)$ at rate $a\theta E\left(t\right)$ with ( $0<\theta <1$ ) and hospitalized persons are assumed to infectious but with a low infectivity (rate ${\eta }_1$ ), so that ( $0<{\eta }_2\le {\eta }_1<1$ ) due to the controlled environment and can recovered at (rate ${\omega }_2$ ) or die from the Covid-19 at (rate

${\delta }_2$ ), the effect of saturated treatment is captured by $g\left(H\left(t\right)\right)=\frac{\alpha H\left(t\right)}{1+\kappa H(\; t\; )}$

where $\alpha$ is the cure rate $\kappa$ is the saturation factor that measures the effect of hospitalized being delayed for treatment. The death compartment $D\left(t\right)$ comprises of persons who died from the disease in both compartment $I\left(t\right)$ and $H\left(t\right)$ over time, with $\mu$ define as the natural death rate in all model compartment, $\sigma$ is the burial rate of decease patient in $H\left(t\right)$. While ${\eta }_2$ is the rate at which dead Covid-19 person infects susceptible individuals, and this connotes the denial of releasing the deceased to family members for burial.

Thus, the force of infection is given as:

$\lambda \left(t\right)=\frac{\gamma \left(I\left(t\right)+{\eta }_{1}H\left(t\right)+{\eta }_{2}D\left(t\right)\right)}{S\left(t\right)+E\left(t\right)+I\left(t\right)+H\left(t\right)+D(\; t\; )}$

$N=S\left(t\right)+E\left(t\right)+I\left(t\right)+H\left(t\right)+R\left(t\right)+D\left(t\right),\forall t\ge 0$ (1)

The SEIHRD epidemic model with saturated treatment function consists of the following equation

$\begin{array}{l}S\left(t\right)=A-\left(\lambda \left(t\right)-\mu \right)S\left(t\right)\hfill \\ E\left(t\right)=\lambda \left(t\right)S\left(t\right)-\left(a+\mu \right)E\left(t\right)\hfill \\ I\left(t\right)=a\left(1-\theta \right)E\left(t\right)-\left({\omega }_1+d_1+\mu \right)I\left(t\right)\hfill \\ H\left(t\right)=a\theta E\left(t\right)-g\left(H\left(t\right)\right)-\left({\omega }_2+d_2+\mu \right)H\left(t\right)\hfill \\ R\left(t\right)=g\left(H\left(t\right)\right)+{\omega }_{1}I\left(t\right)+{\omega }_{2}H\left(t\right)\hfill \\ D\left(t\right)=d_{1}I\left(t\right)+d_{2}H\left(t\right)-\sigma D\left(t\right)\hfill \end{array}$ (2)

where $g\left(H\left(t\right)\right)=\frac{\alpha H\left(t\right)}{1+\kappa H\left(t\right)},\forall \alpha >0,\kappa \ge 0$

The schematic diagram of model system (2) is represented in Figure 1, while the description of state variables and parameter are represented in Table 1.

Figure 1. Flowchart diagram. The rectangular boxes represent the cohort for each compartment and the arrow denotes transition phase for coronavirus model system (3).

Table 1. Description of Variable and Parameter for Covid-19 Model.

However, since the recovered compartment in System (2) does not contribute to Covid-19 disease transmission, $R\left(t\right)$ may be decoupled from other variables in the closed system and the new system becomes:

$\begin{array}{l}S\left(t\right)=A-\lambda \left(t\right)-\mu S\left(t\right)\hfill \\ E\left(t\right)=\lambda \left(t\right)-\left(a+\mu \right)E\left(t\right)\hfill \\ I\left(t\right)=a\left(1-\theta \right)E\left(t\right)-\left({\omega }_1+d_1+\mu \right)I\left(t\right)\hfill \\ H\left(t\right)=a\theta -g\left(H\left(t\right)\right)-\left({\omega }_2+d_2+\mu \right)H\left(t\right)\hfill \\ D\left(t\right)=d_{1}I\left(t\right)+d_{2}H\left(t\right)-\sigma D\left(t\right)\hfill \end{array}$ (3)

$Q\left(t\right)=S\left(t\right)+E\left(t\right)+I\left(t\right)+H(\; t\; )$

3. Model Properties and Analysis

In this section, it’s pertinent to show that all state variables of model (3) are biological meaningful and non-negative for all time (t). In other words, the solution of model (3) with positive initial solution will remain positive for all $t\ge 0$ and this can be established using the following lemma.

Theorem 1. Let the initial solution $S\left(0\right)>0$ , $E\left(0\right)>0$, $I\left(0\right)>0$, $H\left(0\right)>0$, $D\left(0\right)>0$, then the solutions $S\mathrm{,}E\mathrm{,}I\mathrm{,}H\mathrm{,}D$ of model (3) are positive for all $t\ge 0$. Furthermore,

$\underset{t\to \infty }{\lim }\sup Q\left(t\right)\le \frac{A}{\mu }\text{and}\underset{t\to \infty }{\lim }D\left(t\right)\le \frac{A\left(d_1+d_2\right)}{\mu \sigma }$

Lemma 2. The region

$\Omega =\left\{\left(S\left(t\right),E\left(t\right),I\left(t\right),H\left(t\right),D\left(t\right)\right)\in {ℝ}^5:Q\left(t\right)\le \frac{A}{\mu },D\left(t\right)\le \frac{A\left(d_1+d_2\right)}{\mu \sigma }\right\}$

is attracting and positively invariant for model 3.

Proof:

Given that $Q\left(t\right)=S\left(t\right)+E\left(t\right)+I\left(t\right)+H\left(t\right)$ yields

$\frac{\text{d}Q\left(t\right)}{\text{d}t}\le A+\mu Q\left(t\right)$ (4)

Applying the standard comparison theorem, the following solution was obtained:

$Q\left(t\right)\le \frac{A}{\mu }+\left(Q\left(0\right)-\frac{A}{\mu }\right){{\displaystyle \exp }}^{\left(-\mu t\right)},\forall t\ge 0$ (5)

If $Q\left(0\right)\le \frac{A}{\mu }$ then the upper bound of $Q\left(t\right)$ is $\frac{A}{\mu }$ when $t\to \infty$. If $Q\left(0\right)\ge \frac{A}{\mu }$, then $Q\left(t\right)$ decreases to $\frac{A}{\mu }$, when $t\to \infty$ and enters $\Omega$ or approaches asymptotically.

Similarly, since $I\left(t\right)\le Q\left(t\right)$ and $H\left(t\right)\le Q\left(t\right)$ from the last equation of model (3), we have

$D\left(t\right)\le \frac{A\left(d_1+d_2\right)}{\mu }+\sigma D\left(t\right)$ (6)

whose solution is:

$D\left(t\right)=\frac{1}{\sigma }\left(\frac{A\left(d_1+d_2\right)}{\mu }-\sigma \left(\frac{A\left(d_1+d_2\right)}{\mu \sigma }-D\left(0\right)\right){{\displaystyle \exp }}^{-\sigma t}\right)$.

Thus, as $t\to \infty$, $D\left(t\right)\le \frac{A\left(d_1+d_2\right)}{\mu \sigma }$, Hence the feasible solution of model (3)

enters the region $\Omega$. In this region, the model can be considered to be epidemiological and mathematically well posed.

4. Steady State and Stability Analysis

4.1. Steady State

At equilibrium, The disease free equilibrium (E₀) of model (3) is obtained by setting $\left(E\left(t\right)=I\left(t\right)=H\left(t\right)=D\left(t\right)=0\right)$ and the $LHS=0$ which results to:

$E_0=\left(S^{*},E^{*},I^{*},H^{*},D^{*}\right)=\left(\frac{A}{\mu },0,0,0,0\right)$ (7)

4.2. Reproduction Number

The basic reproduction number R₀ also called the control reproduction number of the model (3), measures the average number of new Covid-19 cases generated by a typical infectious individual introduced into a population where a certain fraction is protected. The linear system of Equation (7) is established using the next generation matrix method on model (3) to compute the reproduction number see ( [13] ) with notation F (for the new infection terms) and V (for transfer terms) respectively given by

$F=\left[\begin{array}{cccc}0& \gamma & \gamma {\eta }_1& \gamma {\eta }_2\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],V=\left[\begin{array}{cccc}a+\mu & 0& 0& 0\\ a\left(\theta -1\right)& B_1& 0& 0\\ -a\theta & 0& B_2+\alpha & 0\\ 0& -d_1& -d_2& \sigma \end{array}\right]$

where, $B_1={\omega }_1+d_1+\mu$ and $B_2={\omega }_2+d_2+\mu$.

The reproduction number is defined by the quantity $R_{\text{0}}=\rho \left(F{V}^{-1}\right)$, where $\rho$ is the spectral radius,

$\begin{array}{l}{R}_{\text{0}}=\frac{\gamma a\left[\left(1-\theta \right)\left(B_2+\alpha \right)\left(\sigma +d_1{\eta }_2\right)+\theta B_1\left({\eta }_1\sigma +d_2{\eta }_2\right)\right]}{B_1\sigma \left(B_2+\alpha \right)\left(a+\mu \right)},\\ R_{\text{I}}=\frac{a\gamma \left(1-\theta \right)}{B_1\left(a+\alpha \right)}\mathrm{,}{R}_{\text{H}}=\frac{a\gamma {\eta }_1\theta }{\left(B_2+\alpha \right)\left(a+\mu \right)}\mathrm{,}\\ R_{\text{D}}=\frac{a\gamma {\eta }_2}{\sigma B_1\left(a+\mu \right)}\left(d_1\left(1-\theta \right)+\frac{\theta B_1{d}_2}{B_2+\alpha }\right)\end{array}$ (8)

This implies that by interpretation R₀ is the sum of three constituent reproduction number R_I, R_H and R_D respectively which give the contribution of undetected-infected cases, detected-hospitalized cases and deceased cases to transmission of Covid-19.

4.3. Steady State Analysis

Lemma 3. The Disease Free Equilibrium E₀ of model (3) is locally asymptotically stable if the basic reproduction number $R_0<1$ and unstable if $R_0>1$.

Proof.

To proof lemma (3), the Jacobian of model 3 is obtained as:

$J_0=\left[\begin{array}{ccccc}-\mu & 0& -\gamma & -\gamma {\eta }_1& -\gamma {\eta }_2\\ 0& -a-\mu & \gamma & \gamma {\eta }_1& \gamma {\eta }_2\\ 0& a\left(1-\theta \right)& -\mu -d_1-w_1& 0& 0\\ 0& a\theta & 0& -\mu -d_2-w_2-\alpha & 0\\ 0& 0& d_1& d_2& -\sigma \end{array}\right]$ (9)

Jacobian matrix J₀, have one negative eigenvalues ( $-\mu$ ) and the remaining can be obtained through the characteristics equation below:

$\Gamma \left(\lambda \right)={\lambda }^4+a_1{\lambda }^3+a_2{\lambda }^2+a_3\lambda +a_4$

$\begin{array}{l}{a}_1=\sigma +\mu +a+B_1+B_2\\ a_2=\left[\left(1-R_{\text{H}}\right)\left(B_2+\alpha \right)+\left(1-R_{\text{I}}\right)B_1\right]\left(a+\mu \right)\\ +\left(\alpha +\mu +a+B_1+B_2\right)\sigma +B_1\left(B_2+\alpha \right)\\ a_3=\left[\left(1-R_{\text{I}}\right)B_1\sigma +\left(1-R_{\text{H}}\right)\left(B_1+\sigma \right)\left(B_2+\alpha \right)\right]\left(a+\mu \right)\\ -R_{\text{I}}{B}_1\left(B_2+\alpha +d_1{\eta }_2\right)\left(a+\mu \right)-d_2\gamma \theta \\ a_4=\left(1-R_0\right)\left(B_2+\alpha \right)\left(a+\mu \right)B_1\sigma \end{array}$

From Equation (8) it can be seen that $R_{\text{I}}<R_{\text{0}}$, $R_{\text{H}}<R_{\text{0}}$ and $R_{\text{D}}<R_{\text{0}}$ and this implies that the coefficient $a_2\mathrm{,}{a}_3\mathrm{,}{a}_4$ are positive whenever $R_0<1$, and thus all the coefficients are positive. Further, the criteria of Rough-Hurtwiz for the fourth order polynomial is $a_i>0$ for $i=1,2,3,4$ and $a_1{a}_2{a}_3-a_1^2{a}_4-a_3^2>0$ can be easily satisfied by using the above coefficients. So, the model (3) at the disease free equilibrium is locally asymptotically stable if $R_0<1$.

4.4. Endemic Equilibria

To explore the possibility of existence of the positive equilibria (say, equilibria where at least one of the infected variable is non-negative)

$E_1=\left(S^{\ast },E^{\ast },I^{\ast },H^{\ast },D^{\ast }\right)$ (10)

representing any arbitrary equilibrium of model (3).

By some simple calculation we obtained the following:

$\begin{array}{l}{E}^{\ast }=\frac{\varkappa H^{\ast }}{a\left(H^{\ast }+1\right)\theta }\mathrm{,}{I}^{\ast }=\frac{\varkappa \left(1-\theta \right)H^{\ast }}{B_1\left(\kappa H^{\ast }+1\right)\theta }\\ D^{\ast }=\frac{\left[\varkappa \left(1-\theta \right)d_1+B_1{d}_2\theta \left(\kappa H^{\ast }+1\right)\right]H^{\ast }}{B_1\sigma \left(\kappa H^{\ast }+1\right)\theta }\mathrm{,}\end{array}$

$S^{\ast }=\frac{\left[\left(1-\theta \right)\left(d_1+\sigma \right)\varkappa +B_1\left(\kappa H^{\ast }+1\right)\left(d_2+\sigma \right)\theta a+\sigma B_1\varkappa \right]\left(a+\mu \right)\varkappa H^{\ast }}{{\zeta }_1}$

with

$\begin{array}{l}\varkappa =B_2\kappa H^{\ast }+B_2+\alpha \\ {\zeta }_1=[a\sigma \left(\gamma -B_1-\gamma \theta \right)\varkappa +\gamma B_1\left(\kappa H^{\ast }+1\right)\left({\eta }_1\sigma +d_2\right)\theta a\\ +\gamma B_2{\eta }_2\left(\kappa H^{\ast }+1\right)\left(1-\theta \right)a+\alpha d_1\left(1-\theta \right)a+\sigma \mu \varkappa B_1]a\left(\kappa H^{\ast }+1\right)\theta \end{array}$

and $H^{\ast }$ is the positive solution for the following polynomial:

$a_0{\left(H^{\ast }\right)}^3+a_1{\left(H^{\ast }\right)}^2+a_2{H}^{\ast }+a_3=0$ (11)

with the coefficient of $H^{\ast }$ given as

$\begin{array}{l}{a}_0=\gamma B_2{\kappa }^2\left(a+\mu \right)[B_2\left(1-\theta \right)\left(\gamma \sigma +\mu \sigma +d_1\left({\eta }_2\gamma +\mu \right)\right)\\ +\theta B_1\left(\left(\sigma +d_2\right)\mu +\gamma \left(d_2{\eta }_2+{\eta }_1\sigma \right)\right)-B_1{B}_2\sigma ]\\ a_1\mathrm{<mo>\; =\; </mo>2}{B}_2\kappa \left(B_2+\alpha \right)[\left(1-\theta \right)\left(\sigma +d_1\right)\left(a+\mu \right)\mu +\gamma \left(d_1{\eta }_2+\sigma \right)\left(1-\theta \right)\mu \\ -\left(a+\mu \right)\left(1-R_{\text{0}}\right)B_1\sigma ]+2\kappa B_1\theta \mu B_2[\left(a+\mu \right)\left(d_2+\sigma \right)\\ +\gamma \left({\eta }_2{d}_2+{\eta }_1\sigma \right)]+\kappa B_1\theta \alpha \left(a+\mu \right)\left[\left(d_2{\eta }_1+{\eta }_1\sigma \right)\gamma +\mu \left(d_2+\sigma \right)\right]\\ -A{K}^2{B}_2\theta B_1\left[B_2\left(d_1{\eta }_2+\sigma \right)\left(a+\mu \right)R_{\text{I}}+a\gamma \theta \left(d_2{\eta }_2+{\eta }_1\sigma \right)-B_2\sigma \left(a+\mu \right)\right]\end{array}$

$\begin{array}{l}{a}_2=\left(a+\mu \right){\left(B_2+\alpha \right)}^2\left[\left(d_1+\sigma \right)\left(1-\theta \right)+\gamma \left(d_1{\eta }_2+\sigma \right)\left(1-\theta \right)-B_1\sigma \left(1-R_{\text{0}}\right)\right]\\ -B_1\theta \left[\left(a+\mu \right)\left(\sigma +d_2\right)+\gamma \left(d_2{\eta }_2+\sigma {\eta }_1\right)\left(1+A\kappa a\theta \right)\right]\\ +A\kappa \theta B_1\left[B_2\sigma -\left(a+\mu \right)\sigma \left(B_2+\alpha \right)\left(R_{\text{0}}-1\right)-B_2\left(d_1{\eta }_2+\sigma \right)\left(a+\mu \right)R_{\text{I}}\right]\\ a_3=\gamma A\theta \left(a+\mu \right)\left(B_2+\alpha \right)\left(1-R_{\text{0}}\right)B_1\sigma \end{array}$ (12)

From Equation (12), it follows that a₀ and a₁ are always positive whether $1<R_0>1$ since all parameters are non negative. Thus, the number of possible positive real roots of Equation (11) depends on the signs of a₂ and a₃. This is analyzed using Descartes’ law of signs on polynomial (11) as presented in Table 2 with the last column indicating the numbers of positive root in each scenario.

However, theorem (4) summarizes the conditions for existence of endemic equilibrium points for Covid-19 model (3):

Theorem 4. The Covid-19 model system (3) has:

(i) no endemic equilibrium otherwise, if all the coefficient of Equation (11) are all positive;

(ii) a unique endemic equilibrium, if $R_0>1$;

(iii) a two (2) endemic equilibrium points, if $R_0<1$;

Thus, It is clear from Theorem 2 Case (ii) that the model has a unique endemic equilibrium whenever $R_0>1$. Furthermore, Case (iii) establishes the existence

Table 2. Number of possible positive roots of Equation (11).

of multiple equilibrium points which indicates the possibility of backward bifurcation (where the locally-asymptotically stable DFE co-exists with a locally-asymptotically stable endemic equilibrium (see, for instance, [11] [12] [14] [15] in the model (3) when $R_0<1$. Hence, to find the backward bifurcation, the discriminant of polynomial (11) is set to zero and we have:

$\Delta =\left(-27a_0^2\right)a_3^2+\left(18a_0{a}_1{a}_2-4a_1^3\right)a_3+\left(a_1^2{a}_2^2-4a_0{a}_2^3\right)=0$ (13)

solving for the critical value of R₀ denoted by R_c

$R_{\text{c}}=1-\frac{a_3}{\gamma A\theta \left(a+\mu \right)\left(B_2+\alpha \right)B_1\sigma }$ (14)

and the positive roots of a₃ in polynomial (13) which can be viewed as a quadratic equation which yields

$a_3=\frac{-b-\sqrt{b^2-4ac}}{2a}$ (15)

where $a=-27a_0^2<0$, $b=18a_0{a}_1{a}_2-4a$ and $c=a_1^2{a}_2^2-4a_0{a}_2^3$ when $R_{\text{0}}=R_{\text{c}}$, the largest positive solution of a₃ is substituted into Equation (14) thus,

$\begin{array}{c}{R}_{\text{c}}=1-\frac{a_3}{\gamma A\theta \left(a+\mu \right)\left(B_2+\alpha \right)B_1\sigma }\\ =1-\frac{-18a_0{a}_1{a}_2+4a_1^3-\sqrt{{\left(18a_0{a}_1{a}_2-4a_1^3\right)}^2+108a_0^2\left(a_1^2{a}_2^2-4a_0{a}_2^3\right)}}{-54a_0^2\left(\gamma A\theta \left(a+\mu \right)\left(B_2+\alpha \right)B_1\sigma \right)}\\ =1-\frac{18a_0{a}_1{a}_2-4a_1^3+\sqrt{{\left(18a_0{a}_1{a}_2-4a_1^3\right)}^2+108a_0^2\left(a_1^2{a}_2^2-4a_0{a}_2^3\right)}}{54a_0^2\left(\gamma A\theta \left(a+\mu \right)\left(B_2+\alpha \right)B_1\sigma \right)}\end{array}$ (16)

therefore Equation (16) gives the analytical expression of a threshold value for the reproduction number, and backward bifurcation would occur for values of R_o, such that $R_{\text{c}}<R_0<1$. The associated bifurcation diagram is depicted in Figure 2. thus the following results is established.

Lemma 5. The Covid-19 model system (3) exhibits backward bifurcation whenever case (iii) of theorem 4 holds and that ${\mathcal{R}}_{\text{c}}<{\mathcal{R}}_0<1$.

Proof.

By applying the following simplification and change of variables on model (3). Let $S=x_1$, $E=x_2$, $H=x_3$, $I=x_4$, $D=x_5$ and $\stackrel{\dot{}}{S}=F_1$, $\stackrel{\dot{}}{E}=F_2$, $\stackrel{\dot{}}{H}=F_3$, $\stackrel{\dot{}}{D}=F_5$, $\stackrel{\dot{}}{D}=F_5$ so that $N=x_1+x_2+x_3+x_4+x_5$. Further, by using vector notation $X={\left(x_1\mathrm{,}{x}_2\mathrm{,}{x}_3\mathrm{,}{x}_4\mathrm{,}{x}_5\right)}^{\text{T}}$ model (3) is transformed as follows:

$\begin{array}{l}{F}_1=A-\frac{\gamma \left(x_3+{\eta }_1{x}_4+{\eta }_2{x}_5\right)}{x_1+x_2+x_3+x_4+x_5}-\mu x_1\hfill \\ F_2=\frac{\gamma \left(x_3+{\eta }_1{x}_4+{\eta }_2{x}_5\right)}{x_1+x_2+x_3+x_4+x_5}-\left(a+\mu \right)x_2\hfill \\ F_3=a\left(1-\theta \right)x_2-\left({\omega }_1+d_1+\mu \right)x_3\hfill \\ F_4=a\theta -\frac{\alpha x_4}{1+\kappa x_4}-\left({\omega }_2+d_2+\mu \right)x_4\hfill \\ F_5=d_1{x}_3+d_2{x}_4-\sigma x_5\hfill \end{array}$ (17)

Therefore by considering the case when $R_{\text{o}}=1$ and chosen $\gamma$ as a bifurcation parameter, while solving for ${\gamma }^{\ast }$ from $R_{\text{o}}=1$ gives the critical value as:

${\gamma }^{\ast }=\frac{B_1\sigma \left(B_2+\alpha \right)\left(a+\mu \right)}{a\left(1-\theta \right)\left(B_2+\alpha \right)\left(\sigma +d_1{\eta }_2\right)+a\theta B_1\left({\eta }_1\sigma +d_2{\eta }_2\right)}$ (18)

and evaluating the Jacobian matrix (9)) at E₀ when $\gamma ={\gamma }^{\ast }$ has a simple zero eigenvalue ( ${\lambda }_5=0$ ), whose corresponding right eigenvector $w={\left(w_1\mathrm{,}{w}_2\mathrm{,}{w}_3\mathrm{,}{w}_4\mathrm{,}{w}_5\right)}^{\text{T}}$ and left eigenvector $v={\left(v_1\mathrm{,}{v}_2\mathrm{,}{v}_3\mathrm{,}{v}_4\mathrm{,}{v}_5\right)}^{\text{T}}$ are given by

$\begin{array}{l}{w}_1=-\frac{a+\mu }{\mu },w_2=1,w_3=\frac{a\left(1-\theta \right)}{B_1},w_4=\frac{a\theta }{B_2+\alpha },\\ w_5=\frac{a\theta \left[\left(d_2{\eta }_2+{\eta }_1\sigma \right)-\left(d_1{\eta }_1-d_2\right)\sigma \right]}{\sigma \left(d_1{\eta }_2+\sigma \right)\left(B_2+\alpha \right)}\end{array}$ (19)

furthermore,

$\begin{array}{l}{v}_1=0,v_2=\frac{a\left(1-\theta \right)\left(B_2+\alpha \right)\left(d_1{\eta }_2+\sigma \right)+a{B}_1\theta \left(d_2{\eta }_2+{\eta }_1\sigma \right)}{B_1\left(B_2+\alpha \right)\left(a+\mu \right){\eta }_2},\\ v_3=\frac{\sigma +{\eta }_2{d}_2}{{\eta }_2{B}_1},v_4=\frac{{\eta }_1\sigma +{\eta }_2{d}_2}{{\eta }_2\left(B_2+\alpha \right)},v_5=1\end{array}$

Since v₁ is zero, there is no need to find the derivative of F₁ and the associated non-zero partial derivatives of the right hand side of (17) is given by:

$\begin{array}{l}\frac{{\partial }^2{f}_2}{\partial x_2\partial x_3}=\frac{-\gamma \mu }{A},\frac{{\partial }^2{f}_2}{\partial x_2\partial x_4}=\frac{-\gamma {\eta }_1\mu }{A},\frac{{\partial }^2{f}_2}{\partial x_3^2}=\frac{-2\gamma \mu }{A},\\ \frac{{\partial }^2{f}_2}{\partial x_3\partial x_4}=\frac{-\gamma \mu \left(1+{\eta }_1\right)}{A}\frac{{\partial }^2{f}_2}{\partial x_3\partial x_5}=\frac{-\gamma \mu \left(1+{\eta }_2\right)}{A},\frac{{\partial }^2{f}_2}{\partial x_4^2}=\frac{-2\gamma {\eta }_1\mu }{A},\\ \frac{{\partial }^2{f}_2}{\partial x_4\partial x_5}=\frac{-\gamma \mu \left({\eta }_1+{\eta }_2\right)}{A},\frac{{\partial }^2{f}_2}{\partial x_5^2}=\frac{-2\gamma {\eta }_2\mu }{A},\frac{{\partial }^2{f}_2}{\partial x_3\partial {\gamma }^{\ast }}=1,\\ \frac{{\partial }^2{f}_2}{\partial x_4\partial {\gamma }^{\ast }}={\eta }_1,\frac{{\partial }^2{f}_2}{\partial x_5\partial {\gamma }^{\ast }}={\eta }_2\end{array}$

However, using the express in [16] (in particular Theorem 4), two constants a and b are defined as

$a={\displaystyle \underset{k,i,j=1}{\overset{n}{\sum }}}{v}_k{w}_i{w}_j\frac{{\partial }^2{f}_k}{\partial x_i\partial x_j}\left(0,0\right)b={\displaystyle \underset{k,i=1}{\overset{n}{\sum }}}{v}_k{w}_i\frac{{\partial }^2{f}_k}{\partial x_i\partial {\gamma }^{\ast }}\left(0,0\right)$.

Below are the computation for a and b:

$\begin{array}{l}a={\displaystyle \underset{k,i,j=1}{\overset{5}{\sum }}}{v}_k{w}_i{w}_j\frac{{\partial }^2{f}_k}{\partial x_i\partial x_j}\left(0,0\right)\\ a=-\frac{a\left(1-\theta \right)\mu }{{\eta }_{2}A}\left[\frac{\sigma }{{\beta }_1}+\frac{2a\sigma \left(1-\theta \right)}{{\beta }_1^2}+\frac{a\theta \sigma \left(1+{\eta }_1\right)}{{\beta }_1\left({\beta }_2+\alpha \right)}+\frac{a\theta \left(1+{\eta }_2\right)T}{{\beta }_1}\right]\\ -\frac{a\theta \sigma \mu }{{\eta }_{2}A}\left[\frac{{\eta }_1}{{\beta }_2+\alpha }+\frac{{\eta }_{2}T}{\sigma }+\frac{2a\theta {\eta }_1}{{\left({\beta }_2+\alpha \right)}^2}+\frac{2a\theta {\eta }_2{T}^2}{{\sigma }^2}+\frac{a\theta \left({\eta }_1+{\eta }_2\right)T}{\sigma \left({\beta }_2+\alpha \right)}\right]\\ +\frac{2\left(d_2{\eta }_2+{\eta }_1\sigma \right)a^2{\theta }^2\alpha \kappa }{{\eta }_2{\left({\beta }_2+\alpha \right)}^3}\end{array}$ (20)

$\begin{array}{l}b={\displaystyle \underset{k,i=1}{\overset{5}{\sum }}}{v}_k{w}_i\frac{{\partial }^2{f}_k}{\partial x_i\partial {\gamma }^{\ast }}\left(0,0\right)\\ b=\frac{a\left(1-\theta \right)\left(B_2+\alpha \right)\left(d_1{\eta }_2+\sigma \right)+a{B}_1\theta \left(d_2{\eta }_2+{\eta }_1\sigma \right)}{B_1\left(B_2+\alpha \right)\left(a+\mu \right){\eta }_2}\\ \cdot \left[\frac{a\left(1-\theta \right)}{{\beta }_1}+\frac{a\theta {\eta }_1}{{\beta }_2+\alpha }+\frac{a\theta {\eta }_{2}T}{\sigma }\right]\end{array}$ (21)

where $T=\frac{\left(d_1{\eta }_2+{\eta }_1\sigma \right)+\left(d_2-d_1{\eta }_1\right)\sigma }{\left({\beta }_2+\alpha \right)\left(d_1{\eta }_2+\sigma \right)}$.

From Theorem 4, It is established that the direction of bifurcation is determined by the sign of a and b. Furthermore, since the coefficient of b is always positive ( $b>0$ ), it implies that whenever $a>0$ the model exhibit a backward bifurcation and if $a<0$ the model exhibits a forward bifurcation. The proof is complete $\square$.

Furthermore,

$\begin{array}{l}\frac{2\left(d_2{\eta }_2+{\eta }_1\sigma \right)a^2{\theta }^2\alpha \kappa }{{\eta }_2{\left({\beta }_2+\alpha \right)}^3}\\ >\frac{a\left(1-\theta \right)\mu }{{\eta }_{2}A}[\frac{\sigma }{{\beta }_1}+\frac{2a\sigma \left(1-\theta \right)}{{\beta }_1^2}+\frac{a\theta \sigma \left(1+{\eta }_1\right)}{{\beta }_1\left({\beta }_2+\alpha \right)}+\frac{a\theta \left(1+{\eta }_2\right)T}{{\beta }_1}]\\ +\frac{a\theta \sigma \mu }{{\eta }_{2}A}[\frac{{\eta }_1}{{\beta }_2+\alpha }+\frac{{\eta }_{2}T}{\sigma }+\frac{2a\theta {\eta }_1}{{\left({\beta }_2+\alpha \right)}^2}+\frac{2a\theta {\eta }_2{T}^2}{{\sigma }^2}+\frac{a\theta \left({\eta }_1+{\eta }_2\right)T}{\sigma \left({\beta }_2+\alpha \right)}]\end{array}$ (22)

from Equation (22), we can point out that several parameters when stronger than some level will invoke backward bifurcation. Hence the effect of the infected being delayed for treatment, say $\kappa$, is one of the factors which lead to the backward bifurcation amist ( $a\mathrm{,}\theta$ and $\alpha$ ).

Thus the existence of backward bifurcation is crucial and important in planning on how to control Covid-19 disease, the occurrence of a backward bifurcation at $R_0=1$ makes control more difficult which implies that the disease free equilibrium is not globally stable. One control measure often used is the reduction of susceptibility to infection produced by vaccination along with educational programs such as encouragement of better hygiene and preventive measures [17].

Parameter values use in plotting Figure 3 are: $d_1=10,a=8$ for backward bifurcation and Figure 2 are: $d_1=10,a=0.6$ for Forward Bifurcation. Other parameter values are the same as presented in table [5.2].

Figure 2 and Figure 3 depicts the graphical representation for both forward and backward bifurcations respectively. As the reproduction number R₀ passes through 1, we observe changes in the qualitative behavior of the model system (3). For values of R₀ greater than one, forward bifurcation occurs and this implies that Covid-19 disease will persist in the population. However, the presence of a backward bifurcation exemplify catastrophic effects of Covid-19 disease transmission when (R₀) decreases through 1, the infected population may experience

Figure 2. Backward bifurcation.

Figure 3. Forward bifurcation.

a sudden rise as we have observed in the second wave of Covid-19 cases as already warned by Nigeria minister of health [18]. Also when (R₀) decreases through 1, the level of infection remains high as a result, the standard infection-control measure of lowering (R₀) to below 1 will no longer be viable, since (R₀) needs to be below the critical value ( $R_c=0.9605$ ) to achieve effective Covid-19 infection control to bring about disease eradication. Control measures such as wearing face masks, social distancing, isolation of the infected, and case finding need to be implemented together with vaccination, to bring the population to a disease free state. However, we noted that by reducing the parameter range for backward bifurcation to occur is intrinsic to Covid-19 infection dynamics.

5. Parameter Estimation and Sensitivity Analysis

In this section, we perform parameter estimation and sensitivity analysis to provide some level of reliability to simulation of the proposed model.

5.1. Parameter Estimation

In particular, we first fit the model system (3) with observed Nigeria Covid-19 data from NCDC available at [19] using standard model fitting techniques, such as the least squares regression method to estimate parameters ( $A\mathrm{,}{\eta }_1\mathrm{,}{\eta }_2\mathrm{,}\alpha \mathrm{,}\kappa \mathrm{,}\sigma$ ) without realistic estimates as presented in table [5.2, column 2] alongside parameters with realistic estimates (see [7] [20] [21] ). Furthermore the optimal fit of the model to Covid-19 cumulative data is shown in Figure 4 using the total number of Covid-19 confirm cases in Nigeria from February 28 to December 5, 2020.

5.2. Sensitivity Analysis

Sensitivity analysis of the basic reproduction number of Covid-19 model to each of the parameter value is carried out, with the aim to examine behavioral of the model in response to the changes in parameter. However, contribution of parameters to the basic reproduction number reveals their effectiveness or significance to spread of Covid-19 disease in the populace with respect to time. In the present case, the focus is given to determine how changes in the model parameters impact the effective (R₀). Thus, this is done through the normalized

Figure 4. Data fitting using NCDC data.

forward-sensitive index, and this is defined as:

${\Gamma }_{R_0}^{\aleph }=\frac{\partial R_0}{\partial \Pi }\times \frac{\Pi }{R_0}$ (23)

Using (23), the sensitivity index value calculated to any parameter say ( $\Pi$ ) are shown in the last column of Table 3.

6. Simulations and Results

The sensitivity analysis of the basic reproduction number with respect to the model parameters is carried out in section (5.2) and Covid-19 model has been validated by fitting observed data in section (5.1), Thus we perform the numerical simulation of model system (3) using Runge-Kutta fourth order scheme, to elucidate the qualitative analysis of the model system based on variables and parameter values presented in Table 3.

Using the following initial condition for the state variable $S\left(0\right)=170000000$, $E\left(0\right)=30840297$, $I\left(0\right)=43048$, $H\left(0\right)=21524$, $D\left(0\right)=1485$ for the simulation of our proposed model, and also noting the demographic population of Nigeria to be 201,000,000 people.

In Figure 5 and Figure 7, the evolution of cumulative death are obtained as we perturb the transmission rate of expose to infected individuals (a) and effective contact rate of infectious hospitalizes patients ( ${\eta }_1$ ), while other paramaters remain constant.

At different increasing level of transmission rate from expose to infected individual (a), the cumulative death case depicts in Figure 5 shows the impact of

Table 3. Paramater value and sensitivity index.

more individual becoming symptomatic and this connotes high mortality rate, as a result of certain proportion of individuals $\left(1-\theta \right)$ who self isolate and limited medical resources available for hospitalized patients.

Figure 6 depicts the evolution of hospitalized individual in the population and, this reveals that increase in transmission rate from expose to infected individual (a) implies that, more people becomes infected with the disease leading to increase in number of hospitalized individuals [see $a=0.96$ ] in various isolation centers.

Figure 5. Cumulative death case at a = 0.196, 0.296 and 0.96.

Figure 6. Hospitalized population at a = 0.196, 0.296 and 0.96.

In Figure 7, the impact of increase contact rate within the hospitalized patients is varied ( ${\eta }_1$ ) and this reveals a rapid increase in mortality rate due to saturated treatment or lack of sufficient medical resources available in isolation centers [see ( ${\eta }_1=5$ )], which contribute to increase in cumulative death case.

Also Figure 8 depicts the impact of increasing ( ${\eta }_1$ ) in isolated environment, and this shows an increase in number of hospitalized individuals which can be justified by the fact that increased in saturated treatment level contribute to increase in the number of infected Covid-19 individuals. Having in mind that as result of dissatisfaction and quest for alternative medication some of the patients

Figure 7. Cumulative death case at ${\eta }_1=0.05,1$ and 52.

Figure 8. Hospitalized population at ${\eta }_1=0.05,1$ and 5.

might grew weary and decide to source for other means of survival.

7. Conclusion

In this work, we presented a compartmental, deterministic model to assess the impact of undetected infected individuals and infected hospitalized individuals with saturated treatment on the spread of Covid-19 disease transmission dynamics in Nigeria. The model’s basic properties such as positivity and boundedness of solutions, the basic reproduction number as well as the steady states were determined and analysed. The model reveals continuum of disease-free equilibria which was shown to be locally-asymptotically stable whenever the associated basic reproduction number (denoted by R₀) is less than unity, or globally-asymptotically stable for the value of the basic reproduction number less than the critical value ( $R_c=0.96$ ). The model also reveals that even if the transmission rate for exposed Covid-19 patients is lower, controlling the diseases in Nigeria will prove difficult because of the presence of backward bifurcation. The model fitted well to the NCDC (Nigeria center for disease control) cumulative data from Feb. 28 to Dec. 5, 2020. However, the sensitivity analysis of the model parameters reveals that the most sensitive parameters are the force of infection with susceptibility and the rate of transmission from exposing compartment to the infected. Thus, any control measure that can significantly reduce these parameters will contribute effectively to eradicating the disease from the population. Also, the model simulation reveals that an increase in transmission rate for infected Covid-19 patients yields an increase in cumulative death and rapid spread of the disease in the population. Thus decreasing the saturated treatment alone as shown by decreasing ${\eta }_1$ in Figure 7, Figure 8) is not totally plausible. It is, therefore, necessary to implement other control measures against Covid-19 by adhering strictly to all recommended safety procedures to limit the spread of Covid-19 in Nigeria.

Acknowledgements

Sincere thanks to the members of OALib Journal for their professional performance, and special thanks to managing editor Yaqin YAN for a rare attitude of high quality.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Worldmeters (2020) Reported Cases and Deaths by Country, Territory. https://www.worldometers.info/coronavirus/?utm_campaign=CSauthorbio?#countries
[2]	International Health Regulations (2005) Statement on the 2nd Meeting of the International Health Regulations. https://www.who.int/news/item/30-01-2020
[3]	WHO Health Topics (2020) Coronavirus Symptoms. https://www.who.int/health-topics/coronavirus#tab=tab_3
[4]	Ngonghala, C.N., Iboi, E., Eikenberry, S., Scotch, M., MacIntyre, C.R., Bonds, M.H. and Gumel, A.B. (2020) Mathematical Assessment of the Impact of Non-Pharmaceutical Interventions on Curtailing the 2019 Novel Coronavirus. Mathematical Biosciences, 325, Article ID: 108364. https://doi.org/10.1016/j.mbs.2020.108364
[5]	Center for Disease Control and Prevention (2020) Novel Coronavirus Situation Dashboard. Center for Disease Control and Prevention, Atlanta. https://experience.arcgis.com/experirnce/685d0ace521648f8a5beeeee1b9125cd
[6]	Gumel, A.B., Iboi, E.A., Ngonghala, C.N. and Elbasha E.H. (2020) A Primer on Using Mathematics to Understand COVID-19 Dynamics: Modeling, Analysis and Simulations. Infectious Disease Modelling, 6, 148-168. https://doi.org/10.1016/j.idm.2020.11.005
[7]	Iboi, E., Sharomi, O., Ngonghala, C. and Gumel, A.B. (2020) Mathematical Modelling and Analysis of COVID-19 Pandemic in Nigeria. Mathematical Biosciences and Engineering, 17, 7192-7220. https://doi.org/10.3934/mbe.2020369
[8]	Khan, M.A. and Atangana, A (2020) Modeling the Dynamics of Novel Coronavirus (2019-nCov) with Fractional Derivatives. Alexendria Engineering Journal, 59, 2379-2389. https://doi.org/10.1016/j.aej.2020.02.033
[9]	Lin, Q., Zhao, S., Gao, D., Lou, Y., Yang, S., Salihu, S.M., et al. (2020) A Conceptual Model for the Coronavirus Disease 2019 (COVID-19) Outbreak in Wuhan, China with Individual Reaction and Governmental Action. International Journal of Infectious Diseases, 93, 211-216. https://doi.org/10.1016/j.ijid.2020.02.058
[10]	Ivorra B. and Ferrández M. R. and Vela-Pérez M. and Ramos A. M. (2019) Mathematical Modeling of the Spread of the Coronavirus Disease (COVID-19) Taking into Account the Undetected Infections. The Case of China. Communications in Nonlinear Science and Numerical Simulation, 88, Article No. 105303. https://doi.org/10.1016/j.cnsns.2020.105303
[11]	Zhang, J., Jia, J. and Song, X. (2014) Analysis of an SEIR Model with Saturated Incidence and Saturated Treatment Function. The Scientific World Journal, 2014, Article ID: 910421, 11 p. https://doi.org/10.1155/2014/910421
[12]	Zhang, X. and Liu, X. (2008) Backward Bifurcation of an Epidemic Model with Saturated Treatment Function. Journal of Mathematical Analysis and Application, 348, 433-443. https://doi.org/10.1016/j.jmaa.2008.07.042
[13]	Van Den Driessche, P. and Watmough, J. (2002) Reproduction Numbers and Sub-Threshold Endemic Equilibria for Compartmental Models of Disease Transmission. Mathematical Biosciences, 180, 29-48. https://doi.org/10.1016/S0025-5564(02)00108-6
[14]	Pérez, á.G.C., Avila-Vales, E. and Emilio García-Almeida, G. (2019) Bifurcation Analysis of an SIR Model with Logistic Growth. Nonlinear Incidence and Saturated Treatment, 2019, Article ID: 9876013, 21 p. https://doi.org/10.1155/2019/9876013
[15]	Sylvie, D.D.N. (2019) Modelling the Potential Impact of Limited Hospital Beds on Ebola Virus Disease Dynamics. Ph.D Dissertation, Stellenbosch University, Stellenbosch. https://scholar.sun.ac.za
[16]	Castillo-Chavez, C. and Song, B. (2004), Dynamical Models of Tuberculosis and Their Applications. Mathematical Biosciences and Engineering, 1, 361-404. https://doi.org/10.3934/mbe.2004.1.361
[17]	Fred, B., Castillo-Chavez, C. and Song, B. (2004) Mathematical Models for Communicable Disease (CBMS-NSF Regional Conference Series in Applied Mathematics). Series No. 84, Society for Industrial and Applied Mathematics, Philadelphia, 34-37.
[18]	Reuters (2020) Nigeria Warns of Possible New COVID-19 Wave. https://www.voanews.com/africa/nigeria-warns-possible-new-covid-19-wave
[19]	Nigeria Centre for Disease Control (2021) NCDC Situation Report. Nigeria Centre for Disease Control, Abuja. https://ncdc.gov.ng/diseases/sitreps/?cat=14name=An%20update%20of%20COVID-19%20outbreak%20in%20Nigeria
[20]	Madubueze, C.E., Akabuike, N.M. and Dachollom, S. (2020) The Role of Mathematical Model in Curbing COVID-19 in Nigeria. medRxiv Preprint. https://doi.org/10.1101/2020.07.22.20159210
[21]	Yusuf, T.T. and Idisi, O.I. (2020) Modelling the Transmission Dynamics of HIV and HBV Co-Epidemics: Analysis and Simulation. Mathematical Theory and Modeling, 10, 48-77.