# The COVID-19 Model with Partially Recovered Carriers

Authors: T. S. Faniran; E. A. Bakare; A. O. Falade
Journal: Journal of Applied Mathematics
 (2021)

---

## Abstract

Novel coronavirus (COVID-19) has been spreading and wreaking havoc globally, despite massive efforts by the government and World Health Organization (WHO). Consideration of partially recovered carriers is hypothesized to play a leading role in the persistence of the disease and its introduction to new areas. A model for transmission of COVID-19 by symptomless partially recovered carriers is proposed and analysed. It is shown that key parameters can be identified such that below a threshold level, built on these parameters, the epidemic tends towards extinction, while above another threshold, it tends towards a nontrivial epidemic state. Moreover, optimal control analysis of the model, using Pontryagin’s maximum principle, is performed. The optimal controls are characterized in terms of the optimality system and solved numerically for several scenarios. Numerical simulations and sensitivity analysis of the basic reproduction number,Rc, indicate that the disease is mainly driven by parameters involving the partially recovered carriers rather than symptomatic ones. Moreover, optimal control analysis of the model, using Pontryagin’s maximum principle, is performed. The optimal controls were characterized in terms of the optimality system and solved numerically for several scenarios. Numerical simulations were explored to illustrate our theoretical findings, scenarios were built, and the model predicted that social distancing and treatment of the symptomatic will slow down the epidemic curve and reduce mortality of COVID-19 given that there is an average adherence to social distancing and effective treatment are administered.


---

## Body

## 1. Introduction

Coronavirus is one of the major pathogens that primarily targets the human respiratory system. The previous outbreak of coronaviruses include the severe acute respiratory syndrome (SARS), which occurred in 2003 in Mainland China [1], the Middle East respiratory Syndrome (MERS) in 2012 in Saudi Arabia [2], and the MERS outbreak in 2015 in South Korea [3]. In late December 2019, a cluster of patients was admitted to hospitals with an initial diagnosis of pneumonia of an unknown etiology. These patients were etiologically linked to seafood and wet animal wholesale market in Wuhan, the capital city of the Hubei province, China [4]. Early reports predicted the onset of a potential coronavirus outbreak and estimated the reproduction number to be significantly larger than 1 (ranges from 2.24 to 3.58) [5].A report published in Nature revealed that Chinese health authorities concluded that as of February 7, 2020, there have been 31,161 people who have contracted the infection in China and more than 630 people have died of the infection [6] [http://www.nature.com/articles/d41586-020-00154]. Also, the World Health Organization (WHO) reported 51,174 confirmed cases including 15,384 severe cases and 1,666 death cases in China. At the time of preparing this manuscript, the number of cases (March 7, 2020) has passed 100,000 worldwide, reported by Euronews.The most common symptoms at the onset of COVID-19 illness are fever, cough, and fatigue, while other symptoms include sputum production, headache, hemoptysis, diarrhea, and lymphopenia [7].Mathematical models for dynamics of the transmission and control of COVID-19 have been developed to gain insights on the disease dynamics. Also, as recognized by the WHO [8], mathematical models, especially those who are timely, play a key role in informing evidence-based decisions by health decision and policy makers. Only a few mathematical models have so far been publicly released, to the best of our knowledge. Notably among these studies are [9–18].This study has an important difference from those reported above, in that it presents a mathematical model that considers a proportion of infectious patients who recover partially. These proportions may still be virus carriers. According to a new study developed by Lan et al. [19], they followed four medical professionals ages 30 to 36 years who developed COVID-19 and were treated at Wuhan University’s Zhongnan Hospital in China between Jan. 1 and Feb. 15. All of the individuals recovered, and only one was hospitalized during the illness. The patients were treated with oseltamivir, better known under the brand name of Tamiflu, an antiviral drug. The patients were considered recovered after their symptoms. After recovery, the patients were asked to quarantine themselves at home for five days. They continued to undergo throat swabs for the coronavirus after five days for up to 13 days postrecovery. The results showed that every test between day 5 and day 13 was positive for the virus. Also, a paper published by Pappas [20] gave an instance of a Japanese patient who recovered from COVID-19 and then became ill with the disease for a second time. Pappas [20] explained further that one possibility is that the patient’s immune system did not fight off the virus completely as it began to replicate inside her lung. These findings suggested that at least a proportion of infectious individuals did not recover fully from the disease because a small concentration of the virus might still be in their system which began to replicate. Hence, these groups are classified as partially recovered carriers, as they are virus carriers. This is a new finding that has not been modelled mathematically. Therefore, the present study bridges this knowledge gap.The remaining part of this paper is organized as follows: in Section2, model description and formulation are presented. Detailed stability of disease-free equilibrium points is analysed in Section 3. Sections 4 presents the sensitivity analysis of the model parameters while in Section 5, we provide the optimal control problem, its formulation, optimal strategies, proof of existence, and necessary conditions for optimal control. Section 6 provides the numerical simulations of the formulated model. Finally, Section 7 gives the overall concluding remarks of the study.


## 2. Mathematical Formulation of the COVID-19 Model

In this section, a model for the spread of COVID-19 in human and vector population is formulated. The total human population denoted byNH, is partitioned into five classes, namely, the susceptible individual SH, the exposed individual EH, the infectious individual IH, the partially recovered carriers PC, and the fully recovered individuals RH, so that NH=SH+EH+IH+PC+RH.


### 2.1. Assumptions of the Model

The following assumptions were made in order to formulate the equations of the model:(a)
Containment rate,cm, is introduced into the susceptible population(b)
SARS-COV-2 can persist in the body for at least two weeks after symptoms of the disease clear up [20](c)
Infectious patients whose immune system fights off the virus completely with no single virus left progress to the fully recovered class(d)
Infectious patients whose immune system does not fight off the virus completely (i.e., low level of the virus are still present in their system) progress to the partially recovered carrier class(e)
Reverse transcription polymerase chain reaction (RT-PCR) may not be able to detect the presence of a very low level of the virus in the system of partially recovered carriers at the time that they are discharged from the hospital(f)
Virus begins to replicate inside the lungs of the partially recovered carriers again [20](g)
Partially recovered carriers are contagious but show mild or no symptoms(h)
Individuals in the partially recovered carrier class can also recover from COVID-19, but with a significantly slow rateOur model includes a net inflow of susceptible individuals into a region at a rate,λh per unit time. This parameter includes new births, immigration, and emigration. The susceptible individuals acquire COVID-19 through contact with either the infectious individuals or partially recovered carriers at rates α and β. Also, the susceptible population is reduced by containment rate, cm. The proportion of the infectious individuals become partially recovered carriers at a rate 1−ε. Individuals in the partially recovered carrier class recover fully from the disease at a slower rate ψ2η.Applying the assumptions, nomenclature of parameters, and definitions of variables, the following system of ordinary differential equations is formulated:(1)dSHdt=λh−αSHIH−βSHPC−cmSH,(2)dEHdt=αSHIH+βSHPC−θEH,(3)dIHdt=θEH−γIH−ψ1IH−τIH,(4)dPCdt=1−ετIH−ψ2ηPC,(5)dRHdt=ετIH+ψ1IH+ψ2ηPC,with initial conditions(6)SH0=Sh0>0,EH0=EH0,IH0=IH00>0,PC0=PC00,RH0=RH00>0,where the model parameters are nonnegative.Remark 1.
We restructure and modify the model to have more insightful information. However, it has not in any way changed the research focus.
The first four equations, i.e., (1), (2), (3), and (4), are independent of the compartments RH, i.e., (5). Therefore, after decoupling the equations for RH from models (1), (2), (3), (4), and (5), we devote analyses on the remaining equations of (1), (2), (3), and (4) which becomes(7)dSHdt=λh−αSHIH−βSHPC−cmSH,(8)dEHdt=αSHIH+βSHPC−θEH,(9)dIHdt=θEH−γIH−ψ1IH−τIH,(10)dPCdt=1−ετIH−ψ2ηPC,where the disease-free equilibrium point is E0=SH,EH,IH,PC=λh/cm,0,0,0.


## 2.1. Assumptions of the Model

The following assumptions were made in order to formulate the equations of the model:(a)
Containment rate,cm, is introduced into the susceptible population(b)
SARS-COV-2 can persist in the body for at least two weeks after symptoms of the disease clear up [20](c)
Infectious patients whose immune system fights off the virus completely with no single virus left progress to the fully recovered class(d)
Infectious patients whose immune system does not fight off the virus completely (i.e., low level of the virus are still present in their system) progress to the partially recovered carrier class(e)
Reverse transcription polymerase chain reaction (RT-PCR) may not be able to detect the presence of a very low level of the virus in the system of partially recovered carriers at the time that they are discharged from the hospital(f)
Virus begins to replicate inside the lungs of the partially recovered carriers again [20](g)
Partially recovered carriers are contagious but show mild or no symptoms(h)
Individuals in the partially recovered carrier class can also recover from COVID-19, but with a significantly slow rateOur model includes a net inflow of susceptible individuals into a region at a rate,λh per unit time. This parameter includes new births, immigration, and emigration. The susceptible individuals acquire COVID-19 through contact with either the infectious individuals or partially recovered carriers at rates α and β. Also, the susceptible population is reduced by containment rate, cm. The proportion of the infectious individuals become partially recovered carriers at a rate 1−ε. Individuals in the partially recovered carrier class recover fully from the disease at a slower rate ψ2η.Applying the assumptions, nomenclature of parameters, and definitions of variables, the following system of ordinary differential equations is formulated:(1)dSHdt=λh−αSHIH−βSHPC−cmSH,(2)dEHdt=αSHIH+βSHPC−θEH,(3)dIHdt=θEH−γIH−ψ1IH−τIH,(4)dPCdt=1−ετIH−ψ2ηPC,(5)dRHdt=ετIH+ψ1IH+ψ2ηPC,with initial conditions(6)SH0=Sh0>0,EH0=EH0,IH0=IH00>0,PC0=PC00,RH0=RH00>0,where the model parameters are nonnegative.Remark 1.
We restructure and modify the model to have more insightful information. However, it has not in any way changed the research focus.
The first four equations, i.e., (1), (2), (3), and (4), are independent of the compartments RH, i.e., (5). Therefore, after decoupling the equations for RH from models (1), (2), (3), (4), and (5), we devote analyses on the remaining equations of (1), (2), (3), and (4) which becomes(7)dSHdt=λh−αSHIH−βSHPC−cmSH,(8)dEHdt=αSHIH+βSHPC−θEH,(9)dIHdt=θEH−γIH−ψ1IH−τIH,(10)dPCdt=1−ετIH−ψ2ηPC,where the disease-free equilibrium point is E0=SH,EH,IH,PC=λh/cm,0,0,0.


## 3. Stability Analysis

### 3.1. Local Stability of Disease-Free Equilibrium Solution

Using the next-generation operator approach of Diekmann et al., [21], the effective basic reproduction number associated with disease-free equilibrium, E0, and denoted by Rc, is obtained as(11)Rc=αS0ψ2η+βS01−ετψ2ηBT,where(12)S0=λhcm.Clearly, the disease-free equilibrium is locally asymptotically stable ifRc<1. To see this, we obtain the Jacobian matrix of systems (7), (8), (9), and (10) evaluated at E0:(13)JE0=−cm0−αS0−βS00−θαS0βS00θ−BT000τ−ετ−ψ2η,where(14)BT=γ+ψ1+τ.One of the eigenvalues is−cm. The other three are eigenvalues of the matrix(15)−θαS0βS0θ−BT00τ−ετ−ψ2η,whose characteristic equation is(16)a3x3+a2x2+a1x+a0=0,where(17)a0=ψ2ηBTθ−αS0ψ2ηθ−βθS01−ετ,(18)a1=θψ2η+ψ2ηBT+BTθ−αS0θ,(19)a2=θ+BT+ψ2η,(20)a3=1.The Routh-Hurwitz stability criteria for a3×3 matrix is that all values of the determinant of Hurwitz matrices are positive. The determinants of Hurwitz matrices are:(21)D1=∣a2∣>0,where a2>0 and(22)D2=a2a01a1>0.From (22), we have a1a2−a0>0 or a1a2>a0(23)D3=a2a001a100a2a0.(23) gives a0a1a2−a0>0.Clearly,a2 is positive in (19). Also, a1a2>a0 or a1a2−a0>0, in (17), (18), and (19). Finally, we shall show that a0a1a2−a0>0. Since a1a2−a0 is positive, then, we need to show that a0 is also positive so that a0a1a2−a0>0.(24)a0=ψ2ηBTθ−αS0ψ2ηθ+βθS01−ετ>0,1−αS0ψ2ηθ+βθS01−ετψ2ηBTθ>0,1−Rc>0.Hence,1−Rc>0 if Rc<1. Therefore, a0>0 if Rc<1. Hence, by Routh-Hurwitz theorem [22], all the eigenvalues of the Jacobian matrix JE0 have negative real parts when Rc<1 and the disease-free equilibrium solution is locally asymptotically stable if Rc<1. The following theorem summarizes the above result in what follows:Theorem 2.
The disease-free equilibrium solution is locally asymptotically stable ifRc<1.Remark 3.
The biological meaning of the theorem above is that COVID-19 eradication depends on the initial number of the infectious individuals in the population. This implies that a small invasion of infectious individuals into a completely susceptible population will not lead to an outbreak of the disease.
We present in Figures1, 2, 3, and 4 the time series solution of the susceptible, exposed, infected, and partially recovered population of COVID-19, the prevalence against time for the unforced COVID-19 model with associated wavelet spectrum, and the estimated wavelet spectrum of the first week of the year 2 and year 4.5 for the unforced COVID-19 model.Figure 1
Time series solution of susceptible and exposed population of COVID-19.
(a)(b)Figure 2
Time series solution of infected and partially recovered population of COVID-19.
(a)(b)Figure 3
Prevalence against time for the unforced COVID-19 model with associated wavelet spectrum.Figure 4
The estimated wavelet spectrum at the first week of year 2 and year 4.5 for the unforced COVID-19 model.


### 3.2. Global Stability of Disease-Free Equilibrium Solution

The global asymptotic stability of models (19), (20), (22), and (23) is explored in what follows.Theorem 4.
The disease-free equilibriumE0 is globally asymptotically stable if Rc≤1 and unstable if RC>1.Proof.
Consider the Lyapunov function.(25)L=ψ2ηIH+S0β+S0αψ2η1−ετPC.
Its time derivative is(26)L′=ψ2ηdIHdt+S0β+S0αψ2η1−ετdPCdt,L′=ψ2ηθEH−BTIH+S0β+S0αψ2η1−ετ1−ετIH−ψ2ηPC,L′=S0αψ2ηIH+S0β1−ετIH−ψ2ηBTIH+ψ2ηθEH−S0βψ2ηPC−S0αψ22η2PC1−ετ,=ψ2ηBTIHαS0ψ2η+βS01−ετψ2ηBT−1+ψ2ηθEH−S0βPC−S0αψ2ηPC1−ετ,=ψ2ηBTIHRc−1−ψ2ηS0αψ2ηPC1−ετ+S0βPC−θEH≤ψ2ηBTIHRc−1,≤0ifRc≤1.
Therefore,L′≤0 if RC≤1 with L′=0 if and only if Rc=1 or EH=IH=PC=0. It follows that the largest invariant subset in SH,EH,IH,PC∈T:L′=0 is the singleton E0, and hence by LaSalle’s invariance principle [23], the disease-free equilibrium E0 will be approached by all solution trajectories, and hence, the disease-free equilibrium solution is globally asymptotically stable.Remark 5.
The epidemiological implication of the above theorem is that COVID-19 can be eradicated irrespective of the initial sizes of the subpopulation of the model.


### 3.3. Global Stability of the Endemic Equilibrium Point

The endemic equilibrium solutionE1=SH∗,EH∗,IH∗,PC∗ satisfies the following equations:(27)λh−αSH∗IH∗−βSH∗PC∗−cmSH∗=0,(28)αSH∗IH∗+βSH∗PC∗−θEH∗=0,(29)θEH∗−BTIH∗=0,(30)1−ετIH∗−ψ2ηPC∗=0.FindingSH∗, IH∗, and PC∗ in (27), (28), (29), and (30) gives(31)SH∗=λhαIH∗+βPC∗+cm,(32)IH∗=θEH∗BT,(33)PC∗=1−ετIH∗ψ2η.Adding (27) and (28), we get(34)λh−cmSH∗−θEH∗=0.Substituting (31) and (32) in (33) into (34) gives the following(35)−θcmRcEH∗2λh+cmRC−1EH∗=0.From (35), we can obtain EH∗ to be(36)EH∗=λhθ1−1RC.SubstitutingEH∗ in (31), (32), and (33) gives(37)SH∗=λhcmRC,IH∗=λhBT1−1RC,PC∗=1−ετλhBTψ2η1−1RC.Thus, the endemic equilibrium solution exists wheneverRc>1.Furthermore, we establish the global asymptotic stability of the endemic equilibrium solutions converging to the endemic equilibrium point forRc>1. We shall carry out this, by constructing a suitable Lyapunov function of Goh-Volterra type; see [24]. The result below establishes the global stability of the endemic equilibrium solution E1.Theorem 6.
The unique endemic equilibriumE1 is globally asymptotically stable whenever Rc>1.Proof.
Given the following equations which are satisfied by the endemic equilibrium pointE1,(38)λh=αIH∗+βPC∗SH∗+cmSH∗,(39)αIH∗+βPC∗SH∗=θEH∗,(40)θEH∗=BTIH∗,(41)1−ετ=ψ2γPC∗.
Consider the following Goh-Volterra Lyapunov function(42)V=SH−SH∗−SH∗lnSHSH∗+EH−EH∗−EH∗lnEHEH∗+aIH−IH∗−IH∗lnIHIH∗+bPC−PC∗−PC∗lnPCPC∗,where
a=αSH∗/BT and b=βSH∗/ψ2γ,
with the Lyapunov time derivative obtained as(43)V′=1−SH∗SHSH′+1−EH∗EHEH′+a1−IH∗IHIH′+b1−PC∗PCPC′,V′=1−SH∗SHλh−αIH+βPCSH−cmSH+1−EH∗EHαIH+βPCSH−θEH+a1−IH∗IHθEH−BTIH+b1−PC∗PC1−ετIH−ψ2γPC.
Using (38), we have(44)V′=1−SH∗SHαIH∗+βPC∗SH∗+cmSH∗−αIH+βPCSH−cmSH+1−EH∗EHαIH+βPCSH−θEH+a1−IH∗IHθEH−BTIH+b1−PC∗PC1−ετIH−ψ2γPC.
Further simplification gives(45)V′=αIH∗+βP∗CSH∗+θEH∗+aBTIH∗+bψ2γPC−αSH∗+βPC∗SH∗2SH−αIH+βPCSHEH∗EH−aθEHIH∗IH−b1−ετIHPC∗PC+2cmSH∗−cmSH∗2SH−cmSH+αIH+βPCSH∗−θEH+aθEH−aBTIH+b1−ετIH−bψ2γPC.
Replacinga and b by their values and exploiting (38), (39), (40), and (41) give(46)aθ=αSH∗IH∗EH∗,(47)d1−ετ=βSH∗PC∗IH∗.
Using (38), (39), (40), (41), (46), and (47), we have(48)V′=cmSH∗2−SH∗SH−SHSH∗+3αSH∗IH∗−αIH∗SH∗2SH−αSHIHEH∗EH−αSH∗EHIH∗2EH∗IH+3βSH∗PC∗−βPC∗SH∗2SH−βPCSHEH∗EH−βSH∗IHPC∗2IH∗PC,V′=cmSH∗2−SH∗SH−SHSH∗+αSH∗IH∗3−SH∗SH−SHEH∗IHSH∗EHIH∗−EHIH∗EH∗IH+βSH∗PC∗3−SH∗SH−SHEH∗PCSH∗EHPC∗−IHPC∗IH∗PC.
Using arithmetic-geometric means inequality, i.e.,n−a1+a2+⋯+an≤0, where a1.a2⋯an=1 and a1,a2,⋯,an>0, it follows that V′≤0 with V=0 if and only if SH=SH∗, EH=EH∗, IH=IH∗, and PC=PC∗.
Hence, the largest compact invariant subset of the set whereV′=0 is SH,EH,IH,PC=SH∗,EH∗,IH∗,PC∗ and by classical stability theorem of Lyapunov and LaSalle’s Invariance Principle, it follows that every solution in T approaches E1 for Rc>1 as t→∞.Remark 7.
The epidemiological implication of the above result is that COVID-19 will establish itself wheneverRc>1 in the population.


## 3.1. Local Stability of Disease-Free Equilibrium Solution

Using the next-generation operator approach of Diekmann et al., [21], the effective basic reproduction number associated with disease-free equilibrium, E0, and denoted by Rc, is obtained as(11)Rc=αS0ψ2η+βS01−ετψ2ηBT,where(12)S0=λhcm.Clearly, the disease-free equilibrium is locally asymptotically stable ifRc<1. To see this, we obtain the Jacobian matrix of systems (7), (8), (9), and (10) evaluated at E0:(13)JE0=−cm0−αS0−βS00−θαS0βS00θ−BT000τ−ετ−ψ2η,where(14)BT=γ+ψ1+τ.One of the eigenvalues is−cm. The other three are eigenvalues of the matrix(15)−θαS0βS0θ−BT00τ−ετ−ψ2η,whose characteristic equation is(16)a3x3+a2x2+a1x+a0=0,where(17)a0=ψ2ηBTθ−αS0ψ2ηθ−βθS01−ετ,(18)a1=θψ2η+ψ2ηBT+BTθ−αS0θ,(19)a2=θ+BT+ψ2η,(20)a3=1.The Routh-Hurwitz stability criteria for a3×3 matrix is that all values of the determinant of Hurwitz matrices are positive. The determinants of Hurwitz matrices are:(21)D1=∣a2∣>0,where a2>0 and(22)D2=a2a01a1>0.From (22), we have a1a2−a0>0 or a1a2>a0(23)D3=a2a001a100a2a0.(23) gives a0a1a2−a0>0.Clearly,a2 is positive in (19). Also, a1a2>a0 or a1a2−a0>0, in (17), (18), and (19). Finally, we shall show that a0a1a2−a0>0. Since a1a2−a0 is positive, then, we need to show that a0 is also positive so that a0a1a2−a0>0.(24)a0=ψ2ηBTθ−αS0ψ2ηθ+βθS01−ετ>0,1−αS0ψ2ηθ+βθS01−ετψ2ηBTθ>0,1−Rc>0.Hence,1−Rc>0 if Rc<1. Therefore, a0>0 if Rc<1. Hence, by Routh-Hurwitz theorem [22], all the eigenvalues of the Jacobian matrix JE0 have negative real parts when Rc<1 and the disease-free equilibrium solution is locally asymptotically stable if Rc<1. The following theorem summarizes the above result in what follows:Theorem 2.
The disease-free equilibrium solution is locally asymptotically stable ifRc<1.Remark 3.
The biological meaning of the theorem above is that COVID-19 eradication depends on the initial number of the infectious individuals in the population. This implies that a small invasion of infectious individuals into a completely susceptible population will not lead to an outbreak of the disease.
We present in Figures1, 2, 3, and 4 the time series solution of the susceptible, exposed, infected, and partially recovered population of COVID-19, the prevalence against time for the unforced COVID-19 model with associated wavelet spectrum, and the estimated wavelet spectrum of the first week of the year 2 and year 4.5 for the unforced COVID-19 model.Figure 1
Time series solution of susceptible and exposed population of COVID-19.
(a)(b)Figure 2
Time series solution of infected and partially recovered population of COVID-19.
(a)(b)Figure 3
Prevalence against time for the unforced COVID-19 model with associated wavelet spectrum.Figure 4
The estimated wavelet spectrum at the first week of year 2 and year 4.5 for the unforced COVID-19 model.


## 3.2. Global Stability of Disease-Free Equilibrium Solution

The global asymptotic stability of models (19), (20), (22), and (23) is explored in what follows.Theorem 4.
The disease-free equilibriumE0 is globally asymptotically stable if Rc≤1 and unstable if RC>1.Proof.
Consider the Lyapunov function.(25)L=ψ2ηIH+S0β+S0αψ2η1−ετPC.
Its time derivative is(26)L′=ψ2ηdIHdt+S0β+S0αψ2η1−ετdPCdt,L′=ψ2ηθEH−BTIH+S0β+S0αψ2η1−ετ1−ετIH−ψ2ηPC,L′=S0αψ2ηIH+S0β1−ετIH−ψ2ηBTIH+ψ2ηθEH−S0βψ2ηPC−S0αψ22η2PC1−ετ,=ψ2ηBTIHαS0ψ2η+βS01−ετψ2ηBT−1+ψ2ηθEH−S0βPC−S0αψ2ηPC1−ετ,=ψ2ηBTIHRc−1−ψ2ηS0αψ2ηPC1−ετ+S0βPC−θEH≤ψ2ηBTIHRc−1,≤0ifRc≤1.
Therefore,L′≤0 if RC≤1 with L′=0 if and only if Rc=1 or EH=IH=PC=0. It follows that the largest invariant subset in SH,EH,IH,PC∈T:L′=0 is the singleton E0, and hence by LaSalle’s invariance principle [23], the disease-free equilibrium E0 will be approached by all solution trajectories, and hence, the disease-free equilibrium solution is globally asymptotically stable.Remark 5.
The epidemiological implication of the above theorem is that COVID-19 can be eradicated irrespective of the initial sizes of the subpopulation of the model.


## 3.3. Global Stability of the Endemic Equilibrium Point

The endemic equilibrium solutionE1=SH∗,EH∗,IH∗,PC∗ satisfies the following equations:(27)λh−αSH∗IH∗−βSH∗PC∗−cmSH∗=0,(28)αSH∗IH∗+βSH∗PC∗−θEH∗=0,(29)θEH∗−BTIH∗=0,(30)1−ετIH∗−ψ2ηPC∗=0.FindingSH∗, IH∗, and PC∗ in (27), (28), (29), and (30) gives(31)SH∗=λhαIH∗+βPC∗+cm,(32)IH∗=θEH∗BT,(33)PC∗=1−ετIH∗ψ2η.Adding (27) and (28), we get(34)λh−cmSH∗−θEH∗=0.Substituting (31) and (32) in (33) into (34) gives the following(35)−θcmRcEH∗2λh+cmRC−1EH∗=0.From (35), we can obtain EH∗ to be(36)EH∗=λhθ1−1RC.SubstitutingEH∗ in (31), (32), and (33) gives(37)SH∗=λhcmRC,IH∗=λhBT1−1RC,PC∗=1−ετλhBTψ2η1−1RC.Thus, the endemic equilibrium solution exists wheneverRc>1.Furthermore, we establish the global asymptotic stability of the endemic equilibrium solutions converging to the endemic equilibrium point forRc>1. We shall carry out this, by constructing a suitable Lyapunov function of Goh-Volterra type; see [24]. The result below establishes the global stability of the endemic equilibrium solution E1.Theorem 6.
The unique endemic equilibriumE1 is globally asymptotically stable whenever Rc>1.Proof.
Given the following equations which are satisfied by the endemic equilibrium pointE1,(38)λh=αIH∗+βPC∗SH∗+cmSH∗,(39)αIH∗+βPC∗SH∗=θEH∗,(40)θEH∗=BTIH∗,(41)1−ετ=ψ2γPC∗.
Consider the following Goh-Volterra Lyapunov function(42)V=SH−SH∗−SH∗lnSHSH∗+EH−EH∗−EH∗lnEHEH∗+aIH−IH∗−IH∗lnIHIH∗+bPC−PC∗−PC∗lnPCPC∗,where
a=αSH∗/BT and b=βSH∗/ψ2γ,
with the Lyapunov time derivative obtained as(43)V′=1−SH∗SHSH′+1−EH∗EHEH′+a1−IH∗IHIH′+b1−PC∗PCPC′,V′=1−SH∗SHλh−αIH+βPCSH−cmSH+1−EH∗EHαIH+βPCSH−θEH+a1−IH∗IHθEH−BTIH+b1−PC∗PC1−ετIH−ψ2γPC.
Using (38), we have(44)V′=1−SH∗SHαIH∗+βPC∗SH∗+cmSH∗−αIH+βPCSH−cmSH+1−EH∗EHαIH+βPCSH−θEH+a1−IH∗IHθEH−BTIH+b1−PC∗PC1−ετIH−ψ2γPC.
Further simplification gives(45)V′=αIH∗+βP∗CSH∗+θEH∗+aBTIH∗+bψ2γPC−αSH∗+βPC∗SH∗2SH−αIH+βPCSHEH∗EH−aθEHIH∗IH−b1−ετIHPC∗PC+2cmSH∗−cmSH∗2SH−cmSH+αIH+βPCSH∗−θEH+aθEH−aBTIH+b1−ετIH−bψ2γPC.
Replacinga and b by their values and exploiting (38), (39), (40), and (41) give(46)aθ=αSH∗IH∗EH∗,(47)d1−ετ=βSH∗PC∗IH∗.
Using (38), (39), (40), (41), (46), and (47), we have(48)V′=cmSH∗2−SH∗SH−SHSH∗+3αSH∗IH∗−αIH∗SH∗2SH−αSHIHEH∗EH−αSH∗EHIH∗2EH∗IH+3βSH∗PC∗−βPC∗SH∗2SH−βPCSHEH∗EH−βSH∗IHPC∗2IH∗PC,V′=cmSH∗2−SH∗SH−SHSH∗+αSH∗IH∗3−SH∗SH−SHEH∗IHSH∗EHIH∗−EHIH∗EH∗IH+βSH∗PC∗3−SH∗SH−SHEH∗PCSH∗EHPC∗−IHPC∗IH∗PC.
Using arithmetic-geometric means inequality, i.e.,n−a1+a2+⋯+an≤0, where a1.a2⋯an=1 and a1,a2,⋯,an>0, it follows that V′≤0 with V=0 if and only if SH=SH∗, EH=EH∗, IH=IH∗, and PC=PC∗.
Hence, the largest compact invariant subset of the set whereV′=0 is SH,EH,IH,PC=SH∗,EH∗,IH∗,PC∗ and by classical stability theorem of Lyapunov and LaSalle’s Invariance Principle, it follows that every solution in T approaches E1 for Rc>1 as t→∞.Remark 7.
The epidemiological implication of the above result is that COVID-19 will establish itself wheneverRc>1 in the population.


## 4. Uncertainty and Sensitivity Analysis

### 4.1. Local Sensitivity Analysis

In this section, we carried out sensitivity analysis of parameters of model systems (7), (8), (9), and (10) in order to determine the relative importance of the model parameters on the disease infection. To determine how best to reduce the infection, it is necessary to know the relative importance of the different factors responsible for the infections.Sensitivity indices could be computed numerically so as to figure out parameters that have high impact on basic reproduction numberRc and which of the parameters should be given preferential treatment by intervention strategies.Analytically, sensitivity analysis on all parameters which account for disease dynamics is done using the Chitnis et al. [25] approach; we compute sensitivity indices of the Rc which measures initial disease infection and allows us to measure the relative change in a state variable when a variable changes.The normalized forward sensitivity index of a variable to a parameter is the ratio of the 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. The normalized forward sensitivity index of a variable,u, that depends differentially on a parameter, p, is defined as(49)Npu=∂u∂p×pu,for u≠0.Consequently, we derive analytical expression for the sensitivity index ofRC as(50)NpiRc=∂Rc∂pi×piRc,where pi,i∈ℕ denotes each parameter involved in Rc.Using(51)Rc=αS0ψ2η+βS01−ετψ2ηBT,where S0=λh/cm.(52)BT=γ+ψ1+τ.We compute the sensitivity index of each parameter with respect to theRC, for instance(53)NαRc=∂Rc∂α×αRc=0.999.We have Table1 which summarizes the sensitivity indices of Rc with respect to parameters(54)NλhRc,NcmRc,NγRc,Nψ1Rc,NβRc,NεRc,Nψ2Rc,NτRc.Table 1
Summary of the parameters.
ParameterMeaningValueReferenceαContact rate between susceptible and infectious individuals0.75AssumedβContact rate between susceptible and partially recovered individuals0.75AssumedθTransition rate from exposed to infectious class0.143[15]ψ1Recovery rate of infectious individuals0.33029[15]γCOVID-19-induced death rate1.7826×10−5[15]ψ2Recovery rate of partially recovered carriers[0.000315] (day)−1AssumedτTransition rate from infectious class to partially recovered0.1AssumedεProportion of the infectious individuals who recover fully0.5AssumedϕLoss of immunity rate0.0017AssumedηSlower rate of recovery0.5[0-1]VariedcmContainment rate of susceptible individualsVariable (0-1)AssumedλhNatural birth rate600AssumedContainment rate of susceptible individuals (cm).


#### 4.1.1. Interpretation of Sensitivity Indices Obtained in Table2

The computed sensitivity indices onRc with respect to the involved parameters give insights to the model system proposed. Provided all parameters remain constant, most sensitive parameters are α,λh (the contact rate between susceptible and infectious individuals, natural birth rate). The implication is that the increase in contact rate between susceptible and infectious individuals and natural birth rate have high tendency to increase in COVID-19 epidemic due to positive variation on basic reproduction number Rc. γ (COVID-19-induced rate) is the only parameter that has no effect on the spread of the disease. Interestingly, improvement in slower rate of recovery (η) and increase in transition rate from the infectious disease class to the partially recovered (τ) one are major parameters that reduce the disease epidemic. In the same vein, the sensitivity indexes of cm,ψ1,ψ2 reveal that the increase in containment rate of susceptible individuals, recovery rate of infectious individuals, and recovery rate of partially recovered carriers contribute to the reduction in the spread of COVID-19 infection.Table 2
Numerical values of sensitivity indices ofRc with respect to the parameter involved.
Parameter symbolSensitivity indexα+0.999λh+1.000cm−1.00γ0.00ψ1−0.16β0.049ε0.049η−11709.42ψ2−0.46τ−4652.26In light of this, we recommend that the World health Organization (WHO) should put in place health interventions and strategies that will not only improve the slower recovery rate of the partially recovered individuals but also increase the transition rate from the infectious class to partially recovered one.


### 4.2. Global Sensitivity Analysis

The global sensitivity analysis helps to investigate the change in the output values resulting from changes in all parameter values over the ranges in parameters ([26]. To perform the global sensitivity analysis of the reproduction number Rc and the range of values in Table 3. We aim here to establish the most influential parameters in Rc and also give some insightful ideas from the PRCC plot by running a sample size of 1,000. It was observed that the scatter plot in Figures 5–8 shows a positive relation between Rc and parameters α,β, θ, ψ1, ψ2, τ, ϕ, η, λh, ε, and γ have a negative relation between them and Rc. We understand that the positive relation here implies that a high rate of either of these parameters α,β, θ, ψ1, ψ2, τ, ϕ, and η will generate a higher transmission rate during an outbreak, while the negative relation for λh, ε, and γ will aid in decreasing the severity of the disease-induced death rate. Figure 8(a) presents well defined simulation results of the scatter plots in Figures 5–8 and that of the numerical signs in Table 2. Figure 8(a) also indicates that the most influential parameters that can make the disease spread fast are α,θ, ψ1, and τ. Here, any effort to reduce the relative relevance or impact of α, θ, ψ1, and τ will reduce the spreading rate of COVID-19. Figure 9(b) shows that while Rc (reproduction number) is decreasing with an increase in β, it increases with an increase in ψ2.Table 3
Model parameter descriptions and range.
ParameterMeaningRangeαContact rate between susceptible and infectious individuals0.75–1.5βContact rate between susceptible individuals and partially recovered ones0.75–1.5θTransition rate from the exposed to infectious class0.143–1.0ψ1Recovery rate of infectious individuals0.33029–0.5γCOVID-19-induced death rate0.000017826–0.0010982ψ2Recovery rate of partially recovered carriers0.000315–0.04398τTransition rate from the infectious class to the partially recovered one0.1–0.2384εProportion of the infectious individuals who recover fully0.5–1.0ϕLoss of immunity rate0.0017–1.0cmContainment rate of susceptible individuals(0–1)λhNatural birth rate of human0.038–0.5νSlower rate of recovery0.5–1.5Figure 5
The scatter diagrams for some selected parameters inRC (scatter plots of β, θ, and ψ1, in RC).Figure 6
The scatter diagrams for some selected parameters inRC (scatter plots of ψ2, γ, and ε in RC).Figure 7
The scatter diagrams for some selected parameters inRC (scatter plots of ϕ and ν in RC).Figure 8
The scatter diagrams for some selected parameters inRC (scatter plots of α in RC).Figure 9
Partial rank correlation coefficient (PRCC) of the parameters inRC (Tornado plot of parameters in Rc) in Figure 9(a). The coupling effects of both the recovery rate of partially recovered carriers and the contact rate between susceptible individuals and partially recovered ones on Rc (3D plot of ψ2, β, and Rc)in Figure 9(b).
(a)(b)


## 4.1. Local Sensitivity Analysis

In this section, we carried out sensitivity analysis of parameters of model systems (7), (8), (9), and (10) in order to determine the relative importance of the model parameters on the disease infection. To determine how best to reduce the infection, it is necessary to know the relative importance of the different factors responsible for the infections.Sensitivity indices could be computed numerically so as to figure out parameters that have high impact on basic reproduction numberRc and which of the parameters should be given preferential treatment by intervention strategies.Analytically, sensitivity analysis on all parameters which account for disease dynamics is done using the Chitnis et al. [25] approach; we compute sensitivity indices of the Rc which measures initial disease infection and allows us to measure the relative change in a state variable when a variable changes.The normalized forward sensitivity index of a variable to a parameter is the ratio of the 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. The normalized forward sensitivity index of a variable,u, that depends differentially on a parameter, p, is defined as(49)Npu=∂u∂p×pu,for u≠0.Consequently, we derive analytical expression for the sensitivity index ofRC as(50)NpiRc=∂Rc∂pi×piRc,where pi,i∈ℕ denotes each parameter involved in Rc.Using(51)Rc=αS0ψ2η+βS01−ετψ2ηBT,where S0=λh/cm.(52)BT=γ+ψ1+τ.We compute the sensitivity index of each parameter with respect to theRC, for instance(53)NαRc=∂Rc∂α×αRc=0.999.We have Table1 which summarizes the sensitivity indices of Rc with respect to parameters(54)NλhRc,NcmRc,NγRc,Nψ1Rc,NβRc,NεRc,Nψ2Rc,NτRc.Table 1
Summary of the parameters.
ParameterMeaningValueReferenceαContact rate between susceptible and infectious individuals0.75AssumedβContact rate between susceptible and partially recovered individuals0.75AssumedθTransition rate from exposed to infectious class0.143[15]ψ1Recovery rate of infectious individuals0.33029[15]γCOVID-19-induced death rate1.7826×10−5[15]ψ2Recovery rate of partially recovered carriers[0.000315] (day)−1AssumedτTransition rate from infectious class to partially recovered0.1AssumedεProportion of the infectious individuals who recover fully0.5AssumedϕLoss of immunity rate0.0017AssumedηSlower rate of recovery0.5[0-1]VariedcmContainment rate of susceptible individualsVariable (0-1)AssumedλhNatural birth rate600AssumedContainment rate of susceptible individuals (cm).


### 4.1.1. Interpretation of Sensitivity Indices Obtained in Table2

The computed sensitivity indices onRc with respect to the involved parameters give insights to the model system proposed. Provided all parameters remain constant, most sensitive parameters are α,λh (the contact rate between susceptible and infectious individuals, natural birth rate). The implication is that the increase in contact rate between susceptible and infectious individuals and natural birth rate have high tendency to increase in COVID-19 epidemic due to positive variation on basic reproduction number Rc. γ (COVID-19-induced rate) is the only parameter that has no effect on the spread of the disease. Interestingly, improvement in slower rate of recovery (η) and increase in transition rate from the infectious disease class to the partially recovered (τ) one are major parameters that reduce the disease epidemic. In the same vein, the sensitivity indexes of cm,ψ1,ψ2 reveal that the increase in containment rate of susceptible individuals, recovery rate of infectious individuals, and recovery rate of partially recovered carriers contribute to the reduction in the spread of COVID-19 infection.Table 2
Numerical values of sensitivity indices ofRc with respect to the parameter involved.
Parameter symbolSensitivity indexα+0.999λh+1.000cm−1.00γ0.00ψ1−0.16β0.049ε0.049η−11709.42ψ2−0.46τ−4652.26In light of this, we recommend that the World health Organization (WHO) should put in place health interventions and strategies that will not only improve the slower recovery rate of the partially recovered individuals but also increase the transition rate from the infectious class to partially recovered one.


## 4.1.1. Interpretation of Sensitivity Indices Obtained in Table2

The computed sensitivity indices onRc with respect to the involved parameters give insights to the model system proposed. Provided all parameters remain constant, most sensitive parameters are α,λh (the contact rate between susceptible and infectious individuals, natural birth rate). The implication is that the increase in contact rate between susceptible and infectious individuals and natural birth rate have high tendency to increase in COVID-19 epidemic due to positive variation on basic reproduction number Rc. γ (COVID-19-induced rate) is the only parameter that has no effect on the spread of the disease. Interestingly, improvement in slower rate of recovery (η) and increase in transition rate from the infectious disease class to the partially recovered (τ) one are major parameters that reduce the disease epidemic. In the same vein, the sensitivity indexes of cm,ψ1,ψ2 reveal that the increase in containment rate of susceptible individuals, recovery rate of infectious individuals, and recovery rate of partially recovered carriers contribute to the reduction in the spread of COVID-19 infection.Table 2
Numerical values of sensitivity indices ofRc with respect to the parameter involved.
Parameter symbolSensitivity indexα+0.999λh+1.000cm−1.00γ0.00ψ1−0.16β0.049ε0.049η−11709.42ψ2−0.46τ−4652.26In light of this, we recommend that the World health Organization (WHO) should put in place health interventions and strategies that will not only improve the slower recovery rate of the partially recovered individuals but also increase the transition rate from the infectious class to partially recovered one.


## 4.2. Global Sensitivity Analysis

The global sensitivity analysis helps to investigate the change in the output values resulting from changes in all parameter values over the ranges in parameters ([26]. To perform the global sensitivity analysis of the reproduction number Rc and the range of values in Table 3. We aim here to establish the most influential parameters in Rc and also give some insightful ideas from the PRCC plot by running a sample size of 1,000. It was observed that the scatter plot in Figures 5–8 shows a positive relation between Rc and parameters α,β, θ, ψ1, ψ2, τ, ϕ, η, λh, ε, and γ have a negative relation between them and Rc. We understand that the positive relation here implies that a high rate of either of these parameters α,β, θ, ψ1, ψ2, τ, ϕ, and η will generate a higher transmission rate during an outbreak, while the negative relation for λh, ε, and γ will aid in decreasing the severity of the disease-induced death rate. Figure 8(a) presents well defined simulation results of the scatter plots in Figures 5–8 and that of the numerical signs in Table 2. Figure 8(a) also indicates that the most influential parameters that can make the disease spread fast are α,θ, ψ1, and τ. Here, any effort to reduce the relative relevance or impact of α, θ, ψ1, and τ will reduce the spreading rate of COVID-19. Figure 9(b) shows that while Rc (reproduction number) is decreasing with an increase in β, it increases with an increase in ψ2.Table 3
Model parameter descriptions and range.
ParameterMeaningRangeαContact rate between susceptible and infectious individuals0.75–1.5βContact rate between susceptible individuals and partially recovered ones0.75–1.5θTransition rate from the exposed to infectious class0.143–1.0ψ1Recovery rate of infectious individuals0.33029–0.5γCOVID-19-induced death rate0.000017826–0.0010982ψ2Recovery rate of partially recovered carriers0.000315–0.04398τTransition rate from the infectious class to the partially recovered one0.1–0.2384εProportion of the infectious individuals who recover fully0.5–1.0ϕLoss of immunity rate0.0017–1.0cmContainment rate of susceptible individuals(0–1)λhNatural birth rate of human0.038–0.5νSlower rate of recovery0.5–1.5Figure 5
The scatter diagrams for some selected parameters inRC (scatter plots of β, θ, and ψ1, in RC).Figure 6
The scatter diagrams for some selected parameters inRC (scatter plots of ψ2, γ, and ε in RC).Figure 7
The scatter diagrams for some selected parameters inRC (scatter plots of ϕ and ν in RC).Figure 8
The scatter diagrams for some selected parameters inRC (scatter plots of α in RC).Figure 9
Partial rank correlation coefficient (PRCC) of the parameters inRC (Tornado plot of parameters in Rc) in Figure 9(a). The coupling effects of both the recovery rate of partially recovered carriers and the contact rate between susceptible individuals and partially recovered ones on Rc (3D plot of ψ2, β, and Rc)in Figure 9(b).
(a)(b)


## 5. Analysis of Optimal Control Strategies against COVID-19

There are several measures of prevention and control of infectious diseases, but in this paper, we employ the nonpharmaceutical intervention (NPI), like social distancing,u1, and the nonspecific treatment effort, u2.


### 5.1. Social Distancingu1t

This is a nonpharmaceutical intervention (NPI), also known as physical distancing, which refers to the measures taken to prevent the spread of COVID-19 disease by maintaining a specific physical distance between individuals and reducing the number of time persons come in close contact with one another.


### 5.2. Treatmentu2t

This involves the use of an agent, procedure, or regimen, such as a drug to cure or mitigate the COVID-19 disease. This treatment employed here is a nonspecific regimen.Here, we introduce these time-dependent interventionsu1t,u2t where u1t is the time-dependent preventive effort which we refer to as social distancing, while u2t is the time-dependent treatment effort which is an unspecific regimen to slow the spread of COVID-19. These two control functions are expected to be bounded and Lebesgue integrable on the closed interval 0,T, where T is the duration of time for which we apply our control measures. The COVID-19 models (1), (2), (3), (4), (5), and (7) become(55)dSHdt=λh−α1−u1tSHIH−β1−u1tSHPC−cmSH,(56)dEHdt=α1−u1tSHIH+β1−u1tSHPC−θEH,(57)dIHdt=θEH−γIH−ψ1IH−τIH−u2tIH,(58)dPCdt=1−ετIH−ψ2ηPC,(59)dRHdt=ετIH+ψ1IH+ψ2ηPC+u2tIH,with initial conditions(60)SH0=Sh0>0,EH0=EH0,IH0=IH00>0,PC0=PC00,RH0=RH00>0.The model equation can be written in the vectorial form(61)dSHdt=g1SHt,EHt,IHt,PCt,RHt,(62)dEHdt=g2SHt,EHt,IHt,PCt,RHt,(63)dIHdt=g3SHt,EHt,IHt,PCt,RHt,(64)dPCdt=g4SHt,EHt,IHt,PCt,RHt,(65)dRHdt=g5SHt,EHt,IHt,PCt,RHt.NH0=NH,0>0, for NH=SHt,EHt,IHt,PCt,RHt.dNH/dt=Λh−cmSH−γIH≤Λh−cmNH−γNH where the model parameters are nonnegative. The goal here is to determine the optimal control strategy which minimizes the number of symptomatic infectious humans and the cost of the control measures and to predict the impact of these control measures as a means to advice public health officials on the best control and/or elimination policies. The control functions u1t and u2t are defined on the closed interval 0,T, where 0≤u1t≤1 and 0≤u2t≤1. The objective functional is defined mathematically which is given by(66)Ju1t,u2t=∫0TωIH+c12u12+c22u22dt,whereω, c1, and c2 denote the weight constants (it is a parameter which describes the comparative importance of the two terms in the functional) for the relative costs of IH and the control interventions. Hence, there is need to obtain an optimal control pair u1∗t,u2∗t such that(67)Ju1∗t,u2∗=minu1,u2∈ΩJu1t,u2t,whereΩ=u1t,u2t∣0≤u1t≤1,0≤u2t≤1.


### 5.3. The Optimal Control Problem Statement

Here, we present the optimal control problem for the COVID-19 disease which is stated as follows:(68)Ju1t,u2t=∫0TωIH+c12u12+c22u22dt,subject to the dynamics(69)dSHdt=λh−α1−u1tSHIH−β1−u1tSHPC−cmSH,(70)dEHdt=α1−u1tSHIH+β1−u1tSHPC−θEH,(71)dIHdt=θEH−γIH−ψ1IH−τIH−u2tIH,(72)dPCdt=1−ετIH−ψ2ηPC,(73)dRHdt=ετIH+ψ1IH+ψ2ηPC+u2tIH,with initial conditions(74)SH0=Sh0>0,EH0=EH0,IH0=IH00>0,PC0=PC00,RH0=RH00>0.The terminal or final timeT is fixed and subjected to the control constraints.(75)u1,u2∈Ω,∀t∈0,T.There is a need to solve the following optimal control problem such that we find a controlu1∗t, u2∗t which minimizes the objective functional, that is(76)Ju1∗t,u2∗=minu1,u2∈ΩJu1t,u2t.If there exist suchu1t, u2t, it is referred to as optimal control. The optimal control along side with the corresponding solution gives the optimal control pair(77)SH∗t,EH∗t,IH∗t,PC∗t,RH∗t,u1∗t,u2∗t.The first question that we must address is to confirm whether an optimal control pairSH∗t,EH∗t,IH∗t,PC∗t,RH∗t,u1∗t,u2∗t exists. Thus, this question of existence is settled by the following Lemma.Lemma 8.
(Filippov-Cesari existence theorem). For allt,x∈ℝn+1, define the set(78)Lt,x=fx,u+ζ,gx,u:ζ≤0,u∈U.
Suppose that(1)
Lt,x is convex for every t,x(2)
U is compact(3)
There exists a constantM>0 such that xt≤M for all t∈0,T and all admissible pairs x,u
Then,∃ is an optimal pair x∗t,u∗t where u∗t∈Ω.Lemma 9.
Suppose that there exists a solution of the optimal control problem (68)–(75).Proof.
We rewrite systems (68)–(73) as follows(79)ddtSHEHIHPcRH=g1SHt,EHt,IHt,PCt,RHtg2SHt,EHt,IHt,PCt,RHtg3SHt,EHt,IHt,PCt,RHtg4SHt,EHt,IHt,PCt,RHtg5SHt,EHt,IHt,PCt,RHt=gSHt,EHt,IHt,PCt,RHt,where NH=SHt,EHt,IHt,PCt,RHt. We also define a set Vt,N=fN,ui,t+higN,ui,t:h>0,ui∈U according to the Filippov-Cesari theorem [27], [26].
Here, we first show thatVt,N is a convex for t,N.
Suppose thatz1,z2∈Vt,N, we show that Vt,N is a convex for each t,N by proving that the line connecting z1 and z2 is completely in Vt,N. Thus, we need show thatρzi1+1−ρz2∈Vt,N∀ρ∈0,1. For xi∈Vt,N⇒∃h1,h2≤0 and that the control function u1t,u1t∈U such that zi=fN,ui,t+higN,ui,t:h>0,ui∈U for i=1,2. Let fN,ui,t=ωIH+1/2ciui2, i=1,2,and gN,ui,t=RHSof69−73.
Thus, we have(80)ρfN,u1,u2,t+h1+1−ρfN,u3,u4,t+h2=ρωIH+12c1u12+12c2u22+1−ρωIH+12c1u32+12c2u42+ρh1+1−ρh2=ωIH+ρ12c1u12+12c2u22+1−ρ12c1u32+12c2u42+ρh1+1−ρh2.
Suppose thatu5=ρu12+u22+1−ρu32+u42, it is observed that u5∈U and if we let h3≤0.Y=Therefore, the first part of the convex combination belongs to Vt,N. The second component will be checked as follows:(81)ρgN,u1,u2,t+1−ρgN,u3,u4,t=ρλh−α1−u1tSHIH−β1−u1tSHPC+cmSH+α1−u1tSHIH+β1−u1tSHPC−θEH+θEH−γIH−ψ1IH−τIH−u2tIH+1−ετIH−ψ2ηPC+ετIH+ψ1IH+ψ2ηPC+u2tIH+1−ρλh−α1−u3tSHIH−β1−u3tSHPC=cmSH+α1−u3tSHIH+β1−u3tSHPC−θEH+θEH−γIH−ψ1IH−τIH−u4tIH+1−ετIH−ψ2ηPC+ετIH+ψ1IH+ψ2νPC+u4tIH=Λh−ρα1−u1tSHIH+1−ρα1−u3tSHIH−ρβ1−u1tSHPC+1−ρβ1−u3tSHPC−cmSH+−ρα1−u1tSHIH+1−ρα1−u3tSHIH+ρβ1−u1tSHPC+1−ρβ1−u3t−cmSH−θEH+θEH−γIH−ψ1IH−τIH−ρu2t+1−ρu4tIH+1−ετIH−ψ2ηPC+ετIH+ψ1IH+ψ2ηPC+ρu2t+1−ρu4tIH.
Suppose thatu6=ρu1+u3+1−ρu3+u4. It is observed that u6 belongs to U. Hence, we say finally that the convex combination of ρz1+1−ρz2 is in Vt,N,u1,u2 and U is clearly compact. Secondly, we show that the solution of (69)–(73) is bounded.
For all initial conditions inΔ=SHt,EHt,IHt,PCt,RHt∈ℝ5:SHt+EHt+IHt+PCt+RHt≤1, the trajectory obtained from the initial condition remains in a bounded domain included in Δ.
The total population is given byNH=SHt+EHt+IHt+PCt+RHt with the derivative dNH/dt=ΛhNH−μhNH−γIH.
It is observed that(82)0≤IHt≤NHt,Λh−cmNH−γIH≤dNHdt≤ΛhNH−cmNH−γNH,dNHdt≤ΛhNH−μhNH−γNH.
We have thatNHt≤suptN¯H where NH is the solution of equation(83)dNHdt≤Λh−cmNH−γNH.
Thus,suptNH≤maxNH,0,NH. If NH,0≤NH, then maxtNHt≤NH. If any solution exists, we can find it with the Pontryagin’s maximum principle (PMP). This is done first, by incorporating a time-varying Lagrange multiplier λt, whose elements are referred to as the adjoint or costate variables. We hereby analyse the model equations (68), (69), (70), (71), (72), (73), and (75) with control variables u1t,u2t. We seek an optimal control u1∗t,u2∗ by applying the Pontryagin’s maximum principle (PMP) ([28]) and reformulating (68), (69), (70), (71), (72), (73), and (75) into a problem of minimizing pointwise a Hamiltonian H, with respect to the control variables u1t,u2t . Hence, the Hamiltonian H is defined for all t∈0,T by(84)HSHt,EHt,IHt,PCt,RHt,u1t,u2t,λSHt,λEHt,λIHt,λPCt,λRHt=ωIH+c12u12+u22+λSHλh−α1−u1tSHIH−β1−u1tSHPC−cmSH+λEHα1−u1tSHIH+β1−u1tSHPC−θEH+λIHθEH−γIH−ψ1IH−τIH−u2tIH+λPC1−ετIH−ψ2νPC+λRHετIH+ψ1IH+ψ2ηPC+u2tIH,where λSH, λEH, λIH, λPC, λRH represent the adjoint variables or costate variables. Next, we obtain the following using the Pontryagin maximum principle (PMP) and the existence result of the Filippov-Cesari theorem.Lemma 10.
Suppose an optimal controlu1∗ and u2∗ and the corresponding trajectories SH∗t,EH∗t,IH∗t,PC∗t,andRH∗t of state systems (68), (69), (70), (71), (72), and (73) which minimize Ju1∗,u2∗ over Ω. Then, there exist adjoint variables λSHt, λEHt, λIHt, λPCt, and λRHt which satisfy(85)−dλSHdt=δHδSH=−λH−α1−u1tIH−β1−u1tPC−cm−λEHα1−u1tIH+β1−u1tPC−dλEHdt=δHδEH=−λEH−θ−θλIH−dλIHdt=δHδIH=−ω+λSHα1−u1SHNH−λEHα1−u1SH−λIH−γ−ψ1−τ−u2t−λPC1−ετ−λRHτε+ψ1−u2t−dλPCdt=δHδPC=λSHβ1−u1tSH−λSHβ1−u1tSH−λPC−ψ2η−ψ2ηλPCand with transversality conditions λSHT=λEHT=λIHT=λPCT=λRHT=0 and controls u1∗ and u2∗ satisfy the optimality conditions(86)u1∗=max0,min1,αIH∗SH∗λEH−λSH+βPC∗SH∗λEH−λSHc1,u2∗=max0,min1,λIH−λRHIHc2.Proof.
We obtained the nonlinear system of differential equations for the adjoint variables by differentiating the Hamiltonian functionH, which is evaluated at the optimal control such that the adjoint/costate system is given by(87)−dλSHdt=∂H∂SH=−λH−α1−u1tIH−β1−u1tPC−cm−λEHα1−u1tIH+β1−u1tPC−dλEHdt=∂H∂EH=−λEH−θ−θλIH−dλIHdt=∂H∂IH=−ω+λSHα1−u1SH−λEHα1−u1SH−λIH−γ−ψ1−τ−u2t−λPC1−ετ−λRHτε+ψ1−u2t−dλPCdt=∂H∂PC=λSHβ1−u1tSH−λSHβ1−u1tSH−λPC−ψ2ν−μh−ψ2ηλPC,and with transversality conditions λSHT=λEHT=λIHT=λPCT=λRHT=0 and the controls.
From the optimality condition(88)∂H∂u1=c1u1+λSHαSHIH+βSHPC+λEH−αSHIH−βSHPC=0,we obtain(89)u1∗=αIH∗SH∗λEH−λSH+βPC∗SH∗λEH−λSHc1,∂H∂u2=c2u2∗+λIH−IH∗+λRHIH∗=0.
And we obtain(90)u2∗=λIH−λRHIH∗c2.
Then we have(91)u1∗=max0,min1,αIH∗SH∗λEH−λSH+βPC∗SH∗λEH−λSHNHc1,(92)u2∗=max0,min1,λIH−λRHIH∗c2.
To obtain the optimal control and the corresponding prevalenceIH, we therefore solve the following system(93)dSHdt=λh−α1−u1tSHIH−β1−u1tSHPC−cmSH,dEHdt=α1−u1tSHIH+β1−u1tSHPC−θEH,dIHdt=θEH−γIH−ψ1IH−τIH−u2tIH,dPCdt=1−ετIH−ψ2νPC,dRHdt=ετIH+ψ1IH+ψ2νPC+u2tIH,with initial conditions(94)SH0=Sh0>0,EH0=EH0,IH0=IH00>0,PC0=PC00,RH0=RH00>0,(95)−dλSHdt=g1,(96)−dλEHdt=g2,(97)−dλIHdt=g3,(98)−dλPCdt=g4,(99)−dλRHdt=g5,
where(100)g1=−λH−α1−u1tIH−β1−u1tPC−cm−λEHα1−u1tIH+β1−u1tPC,g2=−λEH−θ−θλIH,g3=−ω+λSHα1−u1SH−λEHα1−u1SH−λIH−γ−ψ1−τ−u2t−λPC1−ετ−λRHτε+ψ1−u2t,g4=λSHβ1−u1tSH−λSHβ1−u1tSH−λPC−ψ2η−ψ2ηλPC,g5=−λRH−ϕ,
and with transversality conditionsλSHT=λEHT=λIHT=λPCT=T=0 where u1∗ and u2∗ are given by (91) and (92). We cannot solve systems (95), (96), (97), (98), and (99) manually but numerical methods must be employed.We hereby present the numerical solution of the optimality system and its corresponding optimal control pairs, the parameter values, and the various scenarios and interpretation in the next section.


## 5.1. Social Distancingu1t

This is a nonpharmaceutical intervention (NPI), also known as physical distancing, which refers to the measures taken to prevent the spread of COVID-19 disease by maintaining a specific physical distance between individuals and reducing the number of time persons come in close contact with one another.


## 5.2. Treatmentu2t

This involves the use of an agent, procedure, or regimen, such as a drug to cure or mitigate the COVID-19 disease. This treatment employed here is a nonspecific regimen.Here, we introduce these time-dependent interventionsu1t,u2t where u1t is the time-dependent preventive effort which we refer to as social distancing, while u2t is the time-dependent treatment effort which is an unspecific regimen to slow the spread of COVID-19. These two control functions are expected to be bounded and Lebesgue integrable on the closed interval 0,T, where T is the duration of time for which we apply our control measures. The COVID-19 models (1), (2), (3), (4), (5), and (7) become(55)dSHdt=λh−α1−u1tSHIH−β1−u1tSHPC−cmSH,(56)dEHdt=α1−u1tSHIH+β1−u1tSHPC−θEH,(57)dIHdt=θEH−γIH−ψ1IH−τIH−u2tIH,(58)dPCdt=1−ετIH−ψ2ηPC,(59)dRHdt=ετIH+ψ1IH+ψ2ηPC+u2tIH,with initial conditions(60)SH0=Sh0>0,EH0=EH0,IH0=IH00>0,PC0=PC00,RH0=RH00>0.The model equation can be written in the vectorial form(61)dSHdt=g1SHt,EHt,IHt,PCt,RHt,(62)dEHdt=g2SHt,EHt,IHt,PCt,RHt,(63)dIHdt=g3SHt,EHt,IHt,PCt,RHt,(64)dPCdt=g4SHt,EHt,IHt,PCt,RHt,(65)dRHdt=g5SHt,EHt,IHt,PCt,RHt.NH0=NH,0>0, for NH=SHt,EHt,IHt,PCt,RHt.dNH/dt=Λh−cmSH−γIH≤Λh−cmNH−γNH where the model parameters are nonnegative. The goal here is to determine the optimal control strategy which minimizes the number of symptomatic infectious humans and the cost of the control measures and to predict the impact of these control measures as a means to advice public health officials on the best control and/or elimination policies. The control functions u1t and u2t are defined on the closed interval 0,T, where 0≤u1t≤1 and 0≤u2t≤1. The objective functional is defined mathematically which is given by(66)Ju1t,u2t=∫0TωIH+c12u12+c22u22dt,whereω, c1, and c2 denote the weight constants (it is a parameter which describes the comparative importance of the two terms in the functional) for the relative costs of IH and the control interventions. Hence, there is need to obtain an optimal control pair u1∗t,u2∗t such that(67)Ju1∗t,u2∗=minu1,u2∈ΩJu1t,u2t,whereΩ=u1t,u2t∣0≤u1t≤1,0≤u2t≤1.


## 5.3. The Optimal Control Problem Statement

Here, we present the optimal control problem for the COVID-19 disease which is stated as follows:(68)Ju1t,u2t=∫0TωIH+c12u12+c22u22dt,subject to the dynamics(69)dSHdt=λh−α1−u1tSHIH−β1−u1tSHPC−cmSH,(70)dEHdt=α1−u1tSHIH+β1−u1tSHPC−θEH,(71)dIHdt=θEH−γIH−ψ1IH−τIH−u2tIH,(72)dPCdt=1−ετIH−ψ2ηPC,(73)dRHdt=ετIH+ψ1IH+ψ2ηPC+u2tIH,with initial conditions(74)SH0=Sh0>0,EH0=EH0,IH0=IH00>0,PC0=PC00,RH0=RH00>0.The terminal or final timeT is fixed and subjected to the control constraints.(75)u1,u2∈Ω,∀t∈0,T.There is a need to solve the following optimal control problem such that we find a controlu1∗t, u2∗t which minimizes the objective functional, that is(76)Ju1∗t,u2∗=minu1,u2∈ΩJu1t,u2t.If there exist suchu1t, u2t, it is referred to as optimal control. The optimal control along side with the corresponding solution gives the optimal control pair(77)SH∗t,EH∗t,IH∗t,PC∗t,RH∗t,u1∗t,u2∗t.The first question that we must address is to confirm whether an optimal control pairSH∗t,EH∗t,IH∗t,PC∗t,RH∗t,u1∗t,u2∗t exists. Thus, this question of existence is settled by the following Lemma.Lemma 8.
(Filippov-Cesari existence theorem). For allt,x∈ℝn+1, define the set(78)Lt,x=fx,u+ζ,gx,u:ζ≤0,u∈U.
Suppose that(1)
Lt,x is convex for every t,x(2)
U is compact(3)
There exists a constantM>0 such that xt≤M for all t∈0,T and all admissible pairs x,u
Then,∃ is an optimal pair x∗t,u∗t where u∗t∈Ω.Lemma 9.
Suppose that there exists a solution of the optimal control problem (68)–(75).Proof.
We rewrite systems (68)–(73) as follows(79)ddtSHEHIHPcRH=g1SHt,EHt,IHt,PCt,RHtg2SHt,EHt,IHt,PCt,RHtg3SHt,EHt,IHt,PCt,RHtg4SHt,EHt,IHt,PCt,RHtg5SHt,EHt,IHt,PCt,RHt=gSHt,EHt,IHt,PCt,RHt,where NH=SHt,EHt,IHt,PCt,RHt. We also define a set Vt,N=fN,ui,t+higN,ui,t:h>0,ui∈U according to the Filippov-Cesari theorem [27], [26].
Here, we first show thatVt,N is a convex for t,N.
Suppose thatz1,z2∈Vt,N, we show that Vt,N is a convex for each t,N by proving that the line connecting z1 and z2 is completely in Vt,N. Thus, we need show thatρzi1+1−ρz2∈Vt,N∀ρ∈0,1. For xi∈Vt,N⇒∃h1,h2≤0 and that the control function u1t,u1t∈U such that zi=fN,ui,t+higN,ui,t:h>0,ui∈U for i=1,2. Let fN,ui,t=ωIH+1/2ciui2, i=1,2,and gN,ui,t=RHSof69−73.
Thus, we have(80)ρfN,u1,u2,t+h1+1−ρfN,u3,u4,t+h2=ρωIH+12c1u12+12c2u22+1−ρωIH+12c1u32+12c2u42+ρh1+1−ρh2=ωIH+ρ12c1u12+12c2u22+1−ρ12c1u32+12c2u42+ρh1+1−ρh2.
Suppose thatu5=ρu12+u22+1−ρu32+u42, it is observed that u5∈U and if we let h3≤0.Y=Therefore, the first part of the convex combination belongs to Vt,N. The second component will be checked as follows:(81)ρgN,u1,u2,t+1−ρgN,u3,u4,t=ρλh−α1−u1tSHIH−β1−u1tSHPC+cmSH+α1−u1tSHIH+β1−u1tSHPC−θEH+θEH−γIH−ψ1IH−τIH−u2tIH+1−ετIH−ψ2ηPC+ετIH+ψ1IH+ψ2ηPC+u2tIH+1−ρλh−α1−u3tSHIH−β1−u3tSHPC=cmSH+α1−u3tSHIH+β1−u3tSHPC−θEH+θEH−γIH−ψ1IH−τIH−u4tIH+1−ετIH−ψ2ηPC+ετIH+ψ1IH+ψ2νPC+u4tIH=Λh−ρα1−u1tSHIH+1−ρα1−u3tSHIH−ρβ1−u1tSHPC+1−ρβ1−u3tSHPC−cmSH+−ρα1−u1tSHIH+1−ρα1−u3tSHIH+ρβ1−u1tSHPC+1−ρβ1−u3t−cmSH−θEH+θEH−γIH−ψ1IH−τIH−ρu2t+1−ρu4tIH+1−ετIH−ψ2ηPC+ετIH+ψ1IH+ψ2ηPC+ρu2t+1−ρu4tIH.
Suppose thatu6=ρu1+u3+1−ρu3+u4. It is observed that u6 belongs to U. Hence, we say finally that the convex combination of ρz1+1−ρz2 is in Vt,N,u1,u2 and U is clearly compact. Secondly, we show that the solution of (69)–(73) is bounded.
For all initial conditions inΔ=SHt,EHt,IHt,PCt,RHt∈ℝ5:SHt+EHt+IHt+PCt+RHt≤1, the trajectory obtained from the initial condition remains in a bounded domain included in Δ.
The total population is given byNH=SHt+EHt+IHt+PCt+RHt with the derivative dNH/dt=ΛhNH−μhNH−γIH.
It is observed that(82)0≤IHt≤NHt,Λh−cmNH−γIH≤dNHdt≤ΛhNH−cmNH−γNH,dNHdt≤ΛhNH−μhNH−γNH.
We have thatNHt≤suptN¯H where NH is the solution of equation(83)dNHdt≤Λh−cmNH−γNH.
Thus,suptNH≤maxNH,0,NH. If NH,0≤NH, then maxtNHt≤NH. If any solution exists, we can find it with the Pontryagin’s maximum principle (PMP). This is done first, by incorporating a time-varying Lagrange multiplier λt, whose elements are referred to as the adjoint or costate variables. We hereby analyse the model equations (68), (69), (70), (71), (72), (73), and (75) with control variables u1t,u2t. We seek an optimal control u1∗t,u2∗ by applying the Pontryagin’s maximum principle (PMP) ([28]) and reformulating (68), (69), (70), (71), (72), (73), and (75) into a problem of minimizing pointwise a Hamiltonian H, with respect to the control variables u1t,u2t . Hence, the Hamiltonian H is defined for all t∈0,T by(84)HSHt,EHt,IHt,PCt,RHt,u1t,u2t,λSHt,λEHt,λIHt,λPCt,λRHt=ωIH+c12u12+u22+λSHλh−α1−u1tSHIH−β1−u1tSHPC−cmSH+λEHα1−u1tSHIH+β1−u1tSHPC−θEH+λIHθEH−γIH−ψ1IH−τIH−u2tIH+λPC1−ετIH−ψ2νPC+λRHετIH+ψ1IH+ψ2ηPC+u2tIH,where λSH, λEH, λIH, λPC, λRH represent the adjoint variables or costate variables. Next, we obtain the following using the Pontryagin maximum principle (PMP) and the existence result of the Filippov-Cesari theorem.Lemma 10.
Suppose an optimal controlu1∗ and u2∗ and the corresponding trajectories SH∗t,EH∗t,IH∗t,PC∗t,andRH∗t of state systems (68), (69), (70), (71), (72), and (73) which minimize Ju1∗,u2∗ over Ω. Then, there exist adjoint variables λSHt, λEHt, λIHt, λPCt, and λRHt which satisfy(85)−dλSHdt=δHδSH=−λH−α1−u1tIH−β1−u1tPC−cm−λEHα1−u1tIH+β1−u1tPC−dλEHdt=δHδEH=−λEH−θ−θλIH−dλIHdt=δHδIH=−ω+λSHα1−u1SHNH−λEHα1−u1SH−λIH−γ−ψ1−τ−u2t−λPC1−ετ−λRHτε+ψ1−u2t−dλPCdt=δHδPC=λSHβ1−u1tSH−λSHβ1−u1tSH−λPC−ψ2η−ψ2ηλPCand with transversality conditions λSHT=λEHT=λIHT=λPCT=λRHT=0 and controls u1∗ and u2∗ satisfy the optimality conditions(86)u1∗=max0,min1,αIH∗SH∗λEH−λSH+βPC∗SH∗λEH−λSHc1,u2∗=max0,min1,λIH−λRHIHc2.Proof.
We obtained the nonlinear system of differential equations for the adjoint variables by differentiating the Hamiltonian functionH, which is evaluated at the optimal control such that the adjoint/costate system is given by(87)−dλSHdt=∂H∂SH=−λH−α1−u1tIH−β1−u1tPC−cm−λEHα1−u1tIH+β1−u1tPC−dλEHdt=∂H∂EH=−λEH−θ−θλIH−dλIHdt=∂H∂IH=−ω+λSHα1−u1SH−λEHα1−u1SH−λIH−γ−ψ1−τ−u2t−λPC1−ετ−λRHτε+ψ1−u2t−dλPCdt=∂H∂PC=λSHβ1−u1tSH−λSHβ1−u1tSH−λPC−ψ2ν−μh−ψ2ηλPC,and with transversality conditions λSHT=λEHT=λIHT=λPCT=λRHT=0 and the controls.
From the optimality condition(88)∂H∂u1=c1u1+λSHαSHIH+βSHPC+λEH−αSHIH−βSHPC=0,we obtain(89)u1∗=αIH∗SH∗λEH−λSH+βPC∗SH∗λEH−λSHc1,∂H∂u2=c2u2∗+λIH−IH∗+λRHIH∗=0.
And we obtain(90)u2∗=λIH−λRHIH∗c2.
Then we have(91)u1∗=max0,min1,αIH∗SH∗λEH−λSH+βPC∗SH∗λEH−λSHNHc1,(92)u2∗=max0,min1,λIH−λRHIH∗c2.
To obtain the optimal control and the corresponding prevalenceIH, we therefore solve the following system(93)dSHdt=λh−α1−u1tSHIH−β1−u1tSHPC−cmSH,dEHdt=α1−u1tSHIH+β1−u1tSHPC−θEH,dIHdt=θEH−γIH−ψ1IH−τIH−u2tIH,dPCdt=1−ετIH−ψ2νPC,dRHdt=ετIH+ψ1IH+ψ2νPC+u2tIH,with initial conditions(94)SH0=Sh0>0,EH0=EH0,IH0=IH00>0,PC0=PC00,RH0=RH00>0,(95)−dλSHdt=g1,(96)−dλEHdt=g2,(97)−dλIHdt=g3,(98)−dλPCdt=g4,(99)−dλRHdt=g5,
where(100)g1=−λH−α1−u1tIH−β1−u1tPC−cm−λEHα1−u1tIH+β1−u1tPC,g2=−λEH−θ−θλIH,g3=−ω+λSHα1−u1SH−λEHα1−u1SH−λIH−γ−ψ1−τ−u2t−λPC1−ετ−λRHτε+ψ1−u2t,g4=λSHβ1−u1tSH−λSHβ1−u1tSH−λPC−ψ2η−ψ2ηλPC,g5=−λRH−ϕ,
and with transversality conditionsλSHT=λEHT=λIHT=λPCT=T=0 where u1∗ and u2∗ are given by (91) and (92). We cannot solve systems (95), (96), (97), (98), and (99) manually but numerical methods must be employed.We hereby present the numerical solution of the optimality system and its corresponding optimal control pairs, the parameter values, and the various scenarios and interpretation in the next section.


## 6. Numerical Results

In this section, we present our numerical findings based on the optimal transmission parameter control measures for the COVID-19 dynamic model. We obtained the optimal control by solving the optimality system which involves the nonlinear differential equations (state equation) and the adjoint or costate equations. We applied the iterative scheme to solve the optimality system, and using the fourth-order Runge-Kutta method, we start to solve the state equations while we guess for the control measures over the simulated time. We solved the adjoint equations by a backward fourth-order Runge-Kutta scheme using the present iterative solutions of the state equation because of the transversality conditions. We employed a convex combination of the previous controls to update the controls and the value from the characterizations of the optimality conditions. We repeat the process, and iterations are truncated if the values of the unknowns at the preceding iteration are near to the present iterations ([29]). A COVID-19 model with preventive and treatment as control measures was investigated to predict the effects of control practices and the transmission of COVID-19. Building various scenarios of the two control measures either one control measure at a time or two control measures at a time, we examine and compare numerical results from simulations using the following scenarios.(1)
Strategy A: employing social distancingu1 without treatment u2=0(2)
Strategy B: treatment of the symptomatic individualsu2 without using social distancing u1=0(3)
Strategy C: employing all the two control measuresu1,u2The parameter values are given in Tables1 and 3, and we obtained the optimal control strategy by the forward-backward sweep Runge-Kutta method of order 4. The weight constants ω=500, c1=200, and c2=200. We ran the simulations for T=1,000 days, and the results from the simulations with parameters in Table 3 are presented in Figures 10, 11, and 12.Figure 10
Simulations showing the effect of social distancing only on infected humanIh in Figure 10(a) and the control profile in Figure 10(b).
(a)(b)Figure 11
Simulations showing the effect of treatment only on infected humanIh in Figure 11(a) and the control profile in Figure 11(b).
(a)(b)Figure 12
Simulations showing the effect of social distancing and treatment on infected humanIh in Figure 12(a) and the control profile in Figure 12(b).
(a)(b)In strategy A, we use only the control measureu1 to optimize the objective functional J, while the control measure on treatment u2 is set to zero. The results in Figure 10(a) reveal that the control strategy provide a decrease in the number of symptomatic human IH as against the increase in the uncontrolled case. The control profile is presented in Figure 10(b), where the optimal social distancing control u1 remains constant till the final time.Applying this strategy B, we employed the treatment control measure onlyu2 while we set the control on social distancing u2 to zero. We observed that Figure 11(a) showed that effective treatment only has significant impact in decreasing COVID-19 incidence in the population, while the control profile in Figure 11(b) showed also that the control remains constant till the final time.In strategy C, we applied the two controlsu1,u2, to optimise the objective functional J. For this strategy, in Figure 12(a), it was observed that the control strategy had significant impact in reducing the number of symptomatic humans as against the increased number of cases under the uncontrolled case. The control profile in Figure 12(b) showed also that the control remains constant till the final time.


## 7. Conclusion

Hence, we presented a COVID-19 model with partially recovered carriers. We investigated the stability of the model which revealed that the disease-free equilibrium is locally and globally asymptotically stable. Scatter plots and the tornado plots of the parameters inRc (reproduction number) revealed that the contact rate between susceptible and infectious individuals (α) has a major impact on the transmission of COVID-19 disease, followed by the transmission rate from the exposed to infectious class, θ. This supports the fact that since the transmission rate is very high, the disease can spread fast in the population. We also observed from our results that the influence of the partially recovered contact rate is lower compared to the contact rate between susceptible and infectious individuals. The scenarios built on the optimal control strategies showed that the use of social distancing and treatment are the best option to control the disease in that they reduce the impact of the epidemic in the community and slows down the epidemic curve. Other interesting modeling works can be found in [30, 31]. The developed models predicted the reduction and control of COVID-19 through incorporating multiple control interventions.



---

*Source: 6406274-2021-07-09.xml*