# A Mathematical Modelling and Analysis of COVID-19 Transmission Dynamics with Optimal Control Strategy

Authors: Abadi Abay Gebremeskel; Hailay Weldegiorgis Berhe; Adugna Temesgen Abay
Journal: Computational and Mathematical Methods in Medicine
 (2022)

---

## Abstract

We proposed a deterministic compartmental model for the transmission dynamics of COVID-19 disease. We performed qualitative and quantitative analysis of the deterministic model concerning the local and global stability of the disease-free and endemic equilibrium points. We found that the disease-free equilibrium is locally asymptotically stable when the basic reproduction number is less than unity, while the endemic equilibrium point becomes locally asymptotically stable if the basic reproduction number is above unity. Furthermore, we derived the global stability of both the disease-free and endemic equilibriums of the system by constructing some Lyapunov functions. IfR0≤1, it is found that the disease-free equilibrium is globally asymptotically stable, while the endemic equilibrium point is globally asymptotically stable when R0>1. The numerical results of the general dynamics are in agreement with the theoretical solutions. We established the optimal control strategy by using Pontryagin’s maximum principle. We performed numerical simulations of the optimal control system to investigate the impact of implementing different combinations of optimal controls in controlling and eradicating COVID-19 disease. From this, a significant difference in the number of cases with and without controls was observed. We observed that the implementation of the combination of the control treatment rate, u2, and the control treatment rate, u3, has shown effective and efficient results in eradicating COVID-19 disease in the community relative to the other strategies.


---

## Body

## 1. Introduction

Since the emergence of the disease in December 2019, the coronavirus disease 2019 has been one of the greatest global threats, causing diverse political, economic, and social crises. It is now officially known as COVID-19 by the World Health Organization (WHO). The disease was announced as a pandemic on March 11, 2020, by WHO [1]. The burden of the disease in developing nations like Africa is very serious. The main causes of transmission of COVID-19 are sneezing, coughing, contacting infected people, touching items used daily, and so on [2].Many countries are involved in taking various measures to prevent the pandemic, such as quarantining of suspected individuals, isolation of the infected individuals, lockdown of the community, contact tracing, mask wearing, and physical distancing in collaboration with the WHO. Even though the trial of the production of the medication was not successful, the production of the vaccination seems effective. Vaccinating the susceptible population has officially started globally, including for developing nations.Mathematical modelling of infectious diseases plays a vital role in understanding the underlying mechanisms by which an infectious disease spreads, in forecasting the long-term effect of an epidemic, and in suggesting intervention methods for controlling an existing disease [2–6]. Currently, several authors have studied a mathematical model of COVID-19 dynamics. The majority of the studies on COVID-19 assumed the proposed model as an autonomous or nonautonomous system. Modelling studies such as by [7–9] assumed a deterministic autonomous dynamical system and others considered fractional-order deterministic system [10–13].There are also more studies that have been carried out on the COVID-19 dynamics system using optimal control strategies [13–19]. In [13], the authors proposed a mathematical model involving lock-down, quarantine, and self-isolating control functions to minimize the burden of the disease on society and the associated costs. The authors in [14] developed a deterministic mathematical model with time variable controls. The results of their study revealed that an optimal implementation of personal protective measures reduced the transmission of COVID-19. Quarantine and treatment of the infected individuals are considered optimal controls in the model stated by Abbasi et al. [15]. In addition, the authors employed an impulsive epidemic model to show the abrupt growth in the population. The model proposed by Kada et al. [16] considered two control strategies in their dynamical system. That is, the treatment of COVID-19-infected individuals and using masks for the sensitive body parts, respectively. They found that implementation of the combined control strategies exhibited an effective strategy in reducing COVID-19 epidemic disease in society compared to implementing the individual strategy.A few others proposed a nonautonomous optimal control system [20–22]. In addition to the prevention and treatment measures being implemented so far, many countries have launched a nationwide vaccination campaign against COVID-19, which is a huge step towards reducing the spread of the disease. Therefore, it is vital to study COVID-19 disease dynamics using mathematical modelling, considering it as an autonomous and nonautonomous system with the implementation of vaccination, treatment, and public health education as constant and continuous controls.Therefore, in this paper, we propose a deterministic compartmental model for the dynamics of COVID-19 in a single host population and investigate the mathematical and numerical analysis of the model. More importantly, it establishes the global stability of the disease-free and endemic equilibria and the implication of vaccination, treatment, and public health education on the disease dynamics. The system is further extended into an optimal control problem with the implementation of continuous controls: vaccination, treatment, and public health education. And finally, the best control measure for controlling the disease is recommended.The next part of the paper is organized as follows. Section2 presents the new dynamics model for COVID-19 disease. In Section 3, we study the qualitative analysis including the basic reproduction number, the existence of the disease-free and endemic equilibria, their local and global stability, and the sensitivity of the parameters. The proposed model is further extended into an optimal control problem and is investigated using optimal control theory in Section 4. The numerical methods employed and the numerical simulation of the system with parameter values are performed in Section 5. Finally, the conclusion and recommendation are briefly presented in Section 6.


## 2. Model Formulation

The human population is divided into four subclasses. These subclasses of human populations are as follows: susceptible populationS, infected population I, infected individuals under intensive care unit (ICU) (Iu), and recovered population R. The total population N at a time t is given by Nt=St+It+Iut+Rt. The susceptible population S is increased by the constant recruitment rate Λ. On the other hand, it is reduced by the rate of transmission β, vaccination rate ν, and natural death rate μ. The infected population I is generated proportional to the susceptible and infected individuals at a rate of β. However, the infected compartment I is decreased by the rate of progress α from I to Iu, rate of recovery γ1, and natural death rate μ. The infected population Iu is increased by the transmission rate α but declined by the recovery rate γ2 and natural death rate μ. Moreover, the recovery compartment R gains population from vaccination of susceptible population S at the rate ν, recovery of infected individual I at the rate γ1, and recovery of infected individual Iu at the rate γ2. The class R declined through a natural death rate μ.The following assumptions are considered in the formulation of the model: the population is homogeneously mixed, the exposed compartment is not considered, the infected populationI is the source of infection, but infected individuals under intensive care unit (IUC) Iu are assumed to be hospitalized and restricted, and people contracted with these groups are assumed to take all forms of protection and are thus not assumed to be the source of infection.Following the population scheme in Figure1, the differential equations that describe the transmission of COVID-19 become
(1)dSdt=Λ−βI+ν+μS,dIdt=βSI−α+γ1+μI,dIudt=αI−γ2+δ+μIu,dRdt=γ1I+γ2Iu+νS−μR,Figure 1
Schematic diagram for COVID-19 model.with initial conditions(2)S0=S0>0,I0=I0≥0,Iu0=I0u≥0,R0=R0≥0.The variableR is not involved in the first three equations of the system Equation (1) and does not influence the dynamics of the first three equations. Hence, we only consider the following reduced system in the analysis that follows:
(3)dSdt=Λ−βI+ν+μS,dIdt=βSI−α+γ1+μI,dIudt=αI−γ2+δ+μIu.Descriptions of the state variables and parameters are illustrated in Abbreviations, and the schematic diagram for the system is illustrated in Figure1.


## 3. Mathematical Analysis of the Model

In this section, we must demonstrate that the system Equation (1) is biologically feasible and realistic whenever all of the system Equation (1) state variables are nonnegative for all time t.Theorem 1.
The solutionS,I, Iu, and R of system Equation (1) with initial conditions S0>0,I0≥0,Iu0≥0, and R0≥0 remain positive for all t>0. Furthermore,
(4)limt⟶∞supN≤Λμ.Proof.
Using a similar approach to [23], we obtain
(5)dSdtS=0=Λ,dIdtI=0=0,dIudtIu=0=αI,dRdtR=0=γ1I+γ2Iu+νS.One can observe that the above rates are all nonnegative on the bounding planeR4+. Then, it is easy to show that the region is positively invariant and attracting [24]. The region of attracting for system Equation (1) can be given as
(6)Π=S,I,Iu,R∈R4+:N≤Λμ.Thus, the closed setΠ is positively invariant and attracting region with respect to system Equation (1). The basic reproduction number, denoted by R0, is the very important quantity for qualitative analysis of the dynamics model. It is stated as the average number of secondary cases produced by a single infected individual introduced into a completely susceptible population [7, 24].Using the next-generation matrix approach, as stated in [25], it is easy to determine the basic reproduction number of system Equation (3). One can easily find the disease-free equilibrium of system Equation (3) as E0=Λ/μ+ν,0,0.Let the state variables of system Equation (3) denoted by x′=S,I,Iut. Thus, system Equation (2) can be rewritten as dxi/dt=Fix−Vix, where Fix is the rate of new infection and Vix is the rate of transfer infection compartment into and out of compartment i, and are, respectively, given by
(7)Fix=βSI0,Vix=α+γ1+μIγ2+δ+μIu−αI.The Jacobian matrices ofFix and Vix at E0 are, respectively,
(8)DFix=F=0βS000,DVix=V=α+γ1+μ0−αγ2+δ+μ.It follows that the spectral radius of the next generation matrixFV− is R0 and given by
(9)R0=αβΛμ+νγ1+μ+αγ2+μ+δ.


### 3.1. Existence of Equilibria and Stability

By straightforward computation, system Equation (3) has the disease-free equilibrium E0=Λ/μ+ν,0,0 and a unique endemic equilibrium, E1=S∗,I∗,Iu∗, where
(10)S∗=α+γ1+μβ,I∗=Λβ−μ+νγ1+μ+αβγ1+μ+α,Iu∗=αΛβ−μ+νγ1+μ+αβμ+νγ1+μ+αγ2+μ+δ.In this case,Λβ>μ+νγ1+μ+α to maintain the feasibility state solutions of system Equation (3).


#### 3.1.1. Local Stability of Disease-Free Equilibrium

Theorem 2.
IfR0<1, then the disease-free equilibrium E0 of system Equation (3) is locally asymptotically stable and is unstable otherwise.Proof.
Linearizing system Equation (3) about E0Λ/μ+ν,0,0, we get the characteristic equation at E0.
(11)−μ+ν+λα+γ1+μ−βΛμ+ν+λγ2+δ+μ+λ=0.
It follows thatλ1=−μ+ν,λ2=α+γ1+μR0−1, and λ3=−γ2+δ+μ. This returns out that the disease-free equilibrium, E0, is locally asymptotically stable.


#### 3.1.2. Global Stability of Disease-Free Equilibrium

Theorem 3.
The disease-free equilibrium,E0, of system Equation (3) is globally asymptotically stable whenever R0≤1 and unstable otherwise.Proof.
We define the following Lyapunov function to investigate global stability ofE0 of system Equation (3).
(12)VI,Iu=A1I+A2Iu.
The time derivatives of Lyapunov is defined by(13)dVdt=A1dIdt+A2dIudt,withtheassuptionA1=A2=1.
SubstitutingdI/dt and dIu/dt from system Equation (3) into Equation (13) gives
(14)dVdt=A1βSI−α+γ1+μI+A2αI−γ2+δ+μIu,≤βS0I−α+γ1+μI+αI−γ2+δ+μIu,=βΛμ+νI−α+γ1+μI−γ2+δ+μIu,=α+γ1+νR0−1I−γ2+δ+μIu.Since all the model parameters and variables are nonnegative, it follows from Equation (14) that V̇≤0 for R0≤1 with V̇=0 holds if and only if I=0 or R0=1, and Iu=0. Therefore, V is Lyapunov function of Π. Thus, the state variables I, Iu⟶0 as t⟶∞. Putting I=Iu=0 in Equation (3) shows that S⟶Λ/μ+ν as t⟶∞. It therefore follows from the LaSalle’s invariance principle [26] that, every solution of system Equation (3) with initial conditions in Π approaches E0 as t⟶∞. Therefore, this completes the proof.


### 3.2. Local Stability of the Endemic Equilibrium

Theorem 4.
The endemic equilibrium,E1, of system Equation (3) is locally asymptotically stable whenever R0>1, and unstable if R0<1.Proof.
The following variational matrix is computed atE1 to determine its local stability. That is,
(15)JE1=−βI∗+ν+μ−βS∗0βI∗βS∗−α+γ1+μ00α−γ2+δ+μ.
Following Equation (16), we found the first eigenvalue λ1=−γ2+δ+μ<0. One can find the other eigenvalues from the following reduced variational matrix:
(16)M=−βI∗+ν+μ−α+γ+μν+μ0.The reduced variational matrix has negative eigenvalues iftraceM<0 and detM>0 [27]. Consequently, we then have
(17)traceM<0,detM>0,and the endemic equilibrium pointE1 is locally asymptotically stable.


#### 3.2.1. Global Stability of Endemic Equilibrium

Theorem 5.
IfR0>1, the endemic equilibrium E1 is locally asymptotically stable in the interior of Π.Proof.
Construct the Lyapunov function,V, in SI-plane as
(18)VS,I=S−S∗−S∗lnSS∗+I−I∗−I∗lnII∗.
At the endemic equilibrium, from system Equation (3), we have
(19)Λ=βI∗S∗+ν+μS∗,βI∗S∗=α+γ1+μI∗.
Computing the time derivative ofVS,I along the solutions of system Equation (3), we obtain
(20)dVdt=1−S∗SdSdt+1−I∗IdIdt=1−S∗SΛ−βI+ν+μS+1−I∗IβIS−α+γ1+μI.
Using Equation (19), we obtain
(21)dVdt=1−S∗SβI∗S∗+ν+μS∗−βI+ν+μS+1−I∗IβIS−α+γ1+μI=−ν+μS−S∗2S+βI∗S∗+βIS∗−βIS−βI∗S∗S∗S+βIS+α+γ1+μI∗−βI∗S−α+γ1+μI=−ν+μS−S∗2S+βI∗S∗−βI∗S∗S∗S+βIS∗−βI∗S−α+γ1+μI+βI∗S∗=−ν+μS−S∗2S+βI∗S∗2−SS∗−S∗S=−ν+μS−S∗2S−βI∗S∗x−2+1x.
Letx=S/S∗. Since the arithmetic mean is greater than or equal to the geometric mean, the function x−2+1/x is nonnegative for all x≥0, and S,I,Iu>0 ensures that V̇≤0. Moreover, the equality V̇=0 holds if and only if S=S∗, and I=I∗. Hence, the largest invariant compact set in S,I,Iu∈Π:V̇=0 is the singleton E1. By the LaSalle’s invariant principle [26], the endemic equilibrium E1, if it exists, is globally asymptotically stable inside Π.


### 3.3. Sensitivity Analysis ofR0

It would be of interest to understand the relative importance of model parameters for the transmission and prevalence of communicable and noncommunicable diseases. This will contribute to identifying the critical model parameters that should be taken into account when considering an intervention strategy. In this section, we use a sensitivity analysis of the basic reproduction number,R0, in relation to the various model parameters to quantify the impact of each parameter on decreasing or increasing R0.Following a similar approach [24, 28], the sensitivity analysis of the basic reproduction number, R0, with respect to the parameter Φ, whose sensitivity is to be determined, is computed by
(22)κΦR0=ΦR0×∂R0∂Φ.The normalized forward sensitivity index ofR0 with respect to the parameters is given by
(23)κΛR0=ΛR0×∂R0∂Λ=1,κβR0=βR0×∂R0∂β=1,κγ1R0=γ1R0×∂R0∂γ1=−γ1α+γ1+μ,κμR0=μR0×∂R0∂μ=−μγ1+μ+αγ2+μ+δ+μ+νγ1+γ2+2μ+α+δμ+νγ1+μ+αγ2+μ+δ,καR0=αR0×∂R0∂α=−αα+γ1+μ.Sensitivity indices of the basic reproduction number,R0, are illustrated in Table 1.Table 1
Sensitivity indices ofR0 calculated at the parameter values presented in Table 2.
ParameterSensitivity indexΛ+1β+1μ-0.3638γ1-0.0080α-0.6359One can observe that the sensitivity indicesμ,α, and γ1 are all negative, while the other sensitivity indices Λ and β remain positive. The reproduction number, R0, is the most sensitive to the parameters Λ,β, and α. The implication of this is that an increase (decrease) in the parameters Λ and β causes a corresponding increase (decrease) in R0. Moreover, an increase (decrease) in the parameter α causes to a corresponding decrease (increase) in R0. On the other hand, our sensitivity analysis shows that R0 is less sensitive to γ1 and μ as presented in Table 1.


## 3.1. Existence of Equilibria and Stability

By straightforward computation, system Equation (3) has the disease-free equilibrium E0=Λ/μ+ν,0,0 and a unique endemic equilibrium, E1=S∗,I∗,Iu∗, where
(10)S∗=α+γ1+μβ,I∗=Λβ−μ+νγ1+μ+αβγ1+μ+α,Iu∗=αΛβ−μ+νγ1+μ+αβμ+νγ1+μ+αγ2+μ+δ.In this case,Λβ>μ+νγ1+μ+α to maintain the feasibility state solutions of system Equation (3).


### 3.1.1. Local Stability of Disease-Free Equilibrium

Theorem 2.
IfR0<1, then the disease-free equilibrium E0 of system Equation (3) is locally asymptotically stable and is unstable otherwise.Proof.
Linearizing system Equation (3) about E0Λ/μ+ν,0,0, we get the characteristic equation at E0.
(11)−μ+ν+λα+γ1+μ−βΛμ+ν+λγ2+δ+μ+λ=0.
It follows thatλ1=−μ+ν,λ2=α+γ1+μR0−1, and λ3=−γ2+δ+μ. This returns out that the disease-free equilibrium, E0, is locally asymptotically stable.


### 3.1.2. Global Stability of Disease-Free Equilibrium

Theorem 3.
The disease-free equilibrium,E0, of system Equation (3) is globally asymptotically stable whenever R0≤1 and unstable otherwise.Proof.
We define the following Lyapunov function to investigate global stability ofE0 of system Equation (3).
(12)VI,Iu=A1I+A2Iu.
The time derivatives of Lyapunov is defined by(13)dVdt=A1dIdt+A2dIudt,withtheassuptionA1=A2=1.
SubstitutingdI/dt and dIu/dt from system Equation (3) into Equation (13) gives
(14)dVdt=A1βSI−α+γ1+μI+A2αI−γ2+δ+μIu,≤βS0I−α+γ1+μI+αI−γ2+δ+μIu,=βΛμ+νI−α+γ1+μI−γ2+δ+μIu,=α+γ1+νR0−1I−γ2+δ+μIu.Since all the model parameters and variables are nonnegative, it follows from Equation (14) that V̇≤0 for R0≤1 with V̇=0 holds if and only if I=0 or R0=1, and Iu=0. Therefore, V is Lyapunov function of Π. Thus, the state variables I, Iu⟶0 as t⟶∞. Putting I=Iu=0 in Equation (3) shows that S⟶Λ/μ+ν as t⟶∞. It therefore follows from the LaSalle’s invariance principle [26] that, every solution of system Equation (3) with initial conditions in Π approaches E0 as t⟶∞. Therefore, this completes the proof.


## 3.1.1. Local Stability of Disease-Free Equilibrium

Theorem 2.
IfR0<1, then the disease-free equilibrium E0 of system Equation (3) is locally asymptotically stable and is unstable otherwise.Proof.
Linearizing system Equation (3) about E0Λ/μ+ν,0,0, we get the characteristic equation at E0.
(11)−μ+ν+λα+γ1+μ−βΛμ+ν+λγ2+δ+μ+λ=0.
It follows thatλ1=−μ+ν,λ2=α+γ1+μR0−1, and λ3=−γ2+δ+μ. This returns out that the disease-free equilibrium, E0, is locally asymptotically stable.


## 3.1.2. Global Stability of Disease-Free Equilibrium

Theorem 3.
The disease-free equilibrium,E0, of system Equation (3) is globally asymptotically stable whenever R0≤1 and unstable otherwise.Proof.
We define the following Lyapunov function to investigate global stability ofE0 of system Equation (3).
(12)VI,Iu=A1I+A2Iu.
The time derivatives of Lyapunov is defined by(13)dVdt=A1dIdt+A2dIudt,withtheassuptionA1=A2=1.
SubstitutingdI/dt and dIu/dt from system Equation (3) into Equation (13) gives
(14)dVdt=A1βSI−α+γ1+μI+A2αI−γ2+δ+μIu,≤βS0I−α+γ1+μI+αI−γ2+δ+μIu,=βΛμ+νI−α+γ1+μI−γ2+δ+μIu,=α+γ1+νR0−1I−γ2+δ+μIu.Since all the model parameters and variables are nonnegative, it follows from Equation (14) that V̇≤0 for R0≤1 with V̇=0 holds if and only if I=0 or R0=1, and Iu=0. Therefore, V is Lyapunov function of Π. Thus, the state variables I, Iu⟶0 as t⟶∞. Putting I=Iu=0 in Equation (3) shows that S⟶Λ/μ+ν as t⟶∞. It therefore follows from the LaSalle’s invariance principle [26] that, every solution of system Equation (3) with initial conditions in Π approaches E0 as t⟶∞. Therefore, this completes the proof.


## 3.2. Local Stability of the Endemic Equilibrium

Theorem 4.
The endemic equilibrium,E1, of system Equation (3) is locally asymptotically stable whenever R0>1, and unstable if R0<1.Proof.
The following variational matrix is computed atE1 to determine its local stability. That is,
(15)JE1=−βI∗+ν+μ−βS∗0βI∗βS∗−α+γ1+μ00α−γ2+δ+μ.
Following Equation (16), we found the first eigenvalue λ1=−γ2+δ+μ<0. One can find the other eigenvalues from the following reduced variational matrix:
(16)M=−βI∗+ν+μ−α+γ+μν+μ0.The reduced variational matrix has negative eigenvalues iftraceM<0 and detM>0 [27]. Consequently, we then have
(17)traceM<0,detM>0,and the endemic equilibrium pointE1 is locally asymptotically stable.


### 3.2.1. Global Stability of Endemic Equilibrium

Theorem 5.
IfR0>1, the endemic equilibrium E1 is locally asymptotically stable in the interior of Π.Proof.
Construct the Lyapunov function,V, in SI-plane as
(18)VS,I=S−S∗−S∗lnSS∗+I−I∗−I∗lnII∗.
At the endemic equilibrium, from system Equation (3), we have
(19)Λ=βI∗S∗+ν+μS∗,βI∗S∗=α+γ1+μI∗.
Computing the time derivative ofVS,I along the solutions of system Equation (3), we obtain
(20)dVdt=1−S∗SdSdt+1−I∗IdIdt=1−S∗SΛ−βI+ν+μS+1−I∗IβIS−α+γ1+μI.
Using Equation (19), we obtain
(21)dVdt=1−S∗SβI∗S∗+ν+μS∗−βI+ν+μS+1−I∗IβIS−α+γ1+μI=−ν+μS−S∗2S+βI∗S∗+βIS∗−βIS−βI∗S∗S∗S+βIS+α+γ1+μI∗−βI∗S−α+γ1+μI=−ν+μS−S∗2S+βI∗S∗−βI∗S∗S∗S+βIS∗−βI∗S−α+γ1+μI+βI∗S∗=−ν+μS−S∗2S+βI∗S∗2−SS∗−S∗S=−ν+μS−S∗2S−βI∗S∗x−2+1x.
Letx=S/S∗. Since the arithmetic mean is greater than or equal to the geometric mean, the function x−2+1/x is nonnegative for all x≥0, and S,I,Iu>0 ensures that V̇≤0. Moreover, the equality V̇=0 holds if and only if S=S∗, and I=I∗. Hence, the largest invariant compact set in S,I,Iu∈Π:V̇=0 is the singleton E1. By the LaSalle’s invariant principle [26], the endemic equilibrium E1, if it exists, is globally asymptotically stable inside Π.


## 3.2.1. Global Stability of Endemic Equilibrium

Theorem 5.
IfR0>1, the endemic equilibrium E1 is locally asymptotically stable in the interior of Π.Proof.
Construct the Lyapunov function,V, in SI-plane as
(18)VS,I=S−S∗−S∗lnSS∗+I−I∗−I∗lnII∗.
At the endemic equilibrium, from system Equation (3), we have
(19)Λ=βI∗S∗+ν+μS∗,βI∗S∗=α+γ1+μI∗.
Computing the time derivative ofVS,I along the solutions of system Equation (3), we obtain
(20)dVdt=1−S∗SdSdt+1−I∗IdIdt=1−S∗SΛ−βI+ν+μS+1−I∗IβIS−α+γ1+μI.
Using Equation (19), we obtain
(21)dVdt=1−S∗SβI∗S∗+ν+μS∗−βI+ν+μS+1−I∗IβIS−α+γ1+μI=−ν+μS−S∗2S+βI∗S∗+βIS∗−βIS−βI∗S∗S∗S+βIS+α+γ1+μI∗−βI∗S−α+γ1+μI=−ν+μS−S∗2S+βI∗S∗−βI∗S∗S∗S+βIS∗−βI∗S−α+γ1+μI+βI∗S∗=−ν+μS−S∗2S+βI∗S∗2−SS∗−S∗S=−ν+μS−S∗2S−βI∗S∗x−2+1x.
Letx=S/S∗. Since the arithmetic mean is greater than or equal to the geometric mean, the function x−2+1/x is nonnegative for all x≥0, and S,I,Iu>0 ensures that V̇≤0. Moreover, the equality V̇=0 holds if and only if S=S∗, and I=I∗. Hence, the largest invariant compact set in S,I,Iu∈Π:V̇=0 is the singleton E1. By the LaSalle’s invariant principle [26], the endemic equilibrium E1, if it exists, is globally asymptotically stable inside Π.


## 3.3. Sensitivity Analysis ofR0

It would be of interest to understand the relative importance of model parameters for the transmission and prevalence of communicable and noncommunicable diseases. This will contribute to identifying the critical model parameters that should be taken into account when considering an intervention strategy. In this section, we use a sensitivity analysis of the basic reproduction number,R0, in relation to the various model parameters to quantify the impact of each parameter on decreasing or increasing R0.Following a similar approach [24, 28], the sensitivity analysis of the basic reproduction number, R0, with respect to the parameter Φ, whose sensitivity is to be determined, is computed by
(22)κΦR0=ΦR0×∂R0∂Φ.The normalized forward sensitivity index ofR0 with respect to the parameters is given by
(23)κΛR0=ΛR0×∂R0∂Λ=1,κβR0=βR0×∂R0∂β=1,κγ1R0=γ1R0×∂R0∂γ1=−γ1α+γ1+μ,κμR0=μR0×∂R0∂μ=−μγ1+μ+αγ2+μ+δ+μ+νγ1+γ2+2μ+α+δμ+νγ1+μ+αγ2+μ+δ,καR0=αR0×∂R0∂α=−αα+γ1+μ.Sensitivity indices of the basic reproduction number,R0, are illustrated in Table 1.Table 1
Sensitivity indices ofR0 calculated at the parameter values presented in Table 2.
ParameterSensitivity indexΛ+1β+1μ-0.3638γ1-0.0080α-0.6359One can observe that the sensitivity indicesμ,α, and γ1 are all negative, while the other sensitivity indices Λ and β remain positive. The reproduction number, R0, is the most sensitive to the parameters Λ,β, and α. The implication of this is that an increase (decrease) in the parameters Λ and β causes a corresponding increase (decrease) in R0. Moreover, an increase (decrease) in the parameter α causes to a corresponding decrease (increase) in R0. On the other hand, our sensitivity analysis shows that R0 is less sensitive to γ1 and μ as presented in Table 1.


## 4. Model with Optimal Control

In this section, we will formulate an optimal control problem and analyze it using four time variable control measures,u1t, u2t, u3t, and u4t, to reduce the disease burden and associated costs. The control u1t represents vaccination of the susceptible population to prevent new infection in society, while the controls u2t and u3t represent efforts to provide appropriate treatments to the infected groups, I and Iu, respectively, to reduce the death rate of individuals due to the disease and the associated treatment costs. And the control u4 represents the public health educational rate. That is, educating the community about the importance of handwashing, social distancing, using face masks, staying at home, etc., through different media outlets, to suppress the spread of COVID-19.In system Equation (1), the vaccination rate ν and the treatment rates γ1 and γ2 were considered as constant parameters. However, in this section, we consider them as time variable and control variables u1t,u2t, and u3t, respectively, as stated above. This is a similar approach in [29, 30]. By imposing the time-dependent control variables in system Equation (1), we formulate the optimal control problem and becomes
(24)dSdt=Λ−βI+u1+u4+μS,dIdt=βSI−α+u2+μI,dIudt=αI−δ+u3+μIu,dRdt=u2I+u3Iu+u1S−μR.We introduce the set of all admissible measurable controls defined by(25)U=u=u1,u2,u3,u4:uiisLebesguemeasurable,0≤ui≤1,t∈0,T,i=1,2,3,4.The timeT is the terminal time. For disease control models, it is important to find the control that minimizes the prevalence and cost of controlling the disease. More specifically, we define the objective functional:
(26)Ju≔∫0Tfxt,utdt=∫0TAI+BIu+a1u122+a2u222+a3u322+a4u422dt,where xt=S,I,Iu,R solves the system Equation (24) for the specified control ut. Moreover, the constants A and B, respectively, represent the positive weights of the infected individuals I and Iu, while the constant ai>0 denotes the weight of costs. Further, a1u12/2 represents the cost of vaccination for the susceptible humans, a1u22/2 denotes the cost of treatment for the infectious human individuals I, a3u32/2 represents the cost of infected people under intensive care unit Iu, and a4u42/2 denotes the cost of preventive measures.The aim of the optimal control problem is determining the controlu∗, such that
(27)Ju∗=minJu:u∈U,subject to Equation (24). If such a control u∗ exists, it is, clearly, known as an optimal control. The optimal control, u∗, associated to the corresponding solution, x∗=S∗,I∗,Iu∗,R∗, gives the optimal control pair x∗,u∗.To achieve this, we need to check the existence of an optimal control for the system Equation (24) and the optimality system.


### 4.1. Existence of the Optimal Control

Theorem 6.
ConsiderJu subject to system Equation (24) with S0,I0,Iu0,R0≥0,0,0,0, then there exists an optimal control u∗ and corresponding x∗ that minimizes Ju.We need to check the following assumptions based on [31] to prove the existence of an optimal control.
(i)
The set of controls and corresponding state variables are nonempty(ii)
The control setU is convex and closed(iii)
The right hand side of the system Equation (24) is bounded by a linear function in the state and controlIntegrand of the system Equation (26) is convex on U and is bounded below by C1u12+u22+u32+u42μ/2−C2, where C1>0,C2>0 and μ>1.Proof.
(i)
The control variableU is nonempty set of measurable function on the finite interval time 0,T in the real numbers(ii)
Clearly, the control setU is convex and closed by definition(iii)
The right hand side of system Equation (24) is, by definition, continuous on the finite time interval 0,T, and can be written as a linear function of the control U with coefficients depending on time and state. In addition, all the state and control variables are bounded on 0,T
To verify this condition, letλ∈0,1,v=v1,v2,v3,v4∈U, and w=w1,w2,w3,w4∈U.(28)fx,1−λv+λw−1−λfx,v+λfx,w=a121−λ2v12+2λ1−λv1w1+λ2w12+a221−λ2v22+2λ1−λv2w2+λ2w22+a321−λ2v32+2λ1−λv3w3+λ2w32+a421−λ2v42+2λ1−λv4w4+λ2w42−1−λa12v12+a22v22+a32v32+a42v42−λa12w12+a22w22+a32w32+a42w42=λ2−λa12v1−w12+a22v2−w22+a32v3−w32+a42v4−w42=aλ2−λ2v−w2≤0.
It follows that(29)fx,1−λv+λw≤1−λfx,v+λfx,w.(30)fx,u=AI+BIu+a1u122+a2u222+a3u322+a4u422≥C1u12+u22+u32+u42μ/2−C2,where C1=Mina1,a2,a3,a4,μ=2 and C2>0. Thus, the condition is justified.


### 4.2. Characterization of the Optimal Controls

The characterization of the optimal control is based on the Pontryagin’s maximum principle [32]. This principle converts the optimal control problem into the problem of minimizing pointwise a Hamiltonian with respect to u. Then, the Hamiltonian (H) is given, for all t∈0,T, by
(31)Hxt,ut,λt=AI+BIu+a1u122+a2u222+a3u322+a4u422+λ1Λ−βI+u1+u4+μS+λ2βSI−α+u2+μI+λ3αI−δ+u3+μIu+λ4u2I+u3Iu+u1S−μR.Following the approach [30], if the control u∗ and corresponding state x∗ are the optimal control pair, then necessarily there exists a nonzero adjoint vector function λt fulfilling the following equality.
(32)dxtdt=∂Hxt,ut,λt∂λi,controlsystem,(33)dλtdt=−∂Hxt,ut,λt∂xi,adjointcondition,(34)∂Hxt,ut,λt∂ui=0,optimalitycondition,(35)λT=0,transversalitycondition.Finally, it follows that(36)u∗=0,if∂H∂ui<0,0≤ui∗≤1,if∂H∂ui=0,1,if∂H∂ui>0.Theorem 7 (see [29–31]).
Consider the optimal controlu∗t and corresponding trajectory x∗t with t∈0,T for the system Equation (24), then it is necessary that there exist nontrivial adjoint vector function λt satisfying
(37)dλ1dt=βI+u1+u4+μλ1−βIλ2−u1λ4,dλ2dt=βSλ1−βS−α+μ+u2λ2−αλ3−u2λ4−A,dλ3dt=δ+μ+u3λ3−u3λ4−B,dλ4dt=μλ4,λi=0,i=1,2,3,4.
Hence,(38)Hxt,u∗t,λt=minHxt,ut,λt:u∈U.
Furthermore, the compact form of the optimal control,u∗, is given by
(39)u1∗=max0,min1,λ1−λ4S∗a1,u2∗=max0,min1,λ2−λ4I∗a2,(40)u3∗=max0,min1,λ3−λ4Iu∗a3,u4∗=max0,min1,λ1S∗a4.Proof.
We find the adjoint system Equation (37) by differentiating the Hamiltonian (H) in Equation (31) with respect to x=S,I,Iu,R. More precisely, dλt/dt=−∂Hxt,ut,λt/∂xi. We assume that the final state xT is free. Hence, the transversality condition holds that λi=0,i=1,2,3,4. Moreover, the optimal controls, u1∗,u2∗,u3∗,andu4∗, are obtained based on the optimality condition ∂Hx,u,λ/∂ui=0.
Thus, the controls are become as follows:(41)u1∗=0,ifλ1−λ4S∗a1<0,u1∗,if0≤λ1−λ4S∗a1≤1,1,ifλ1−λ4S∗a1>1,u2∗=0,ifλ2−λ4I∗a2<0,u2∗,if0≤λ2−λ4I∗a2≤1,1,ifλ2−λ4I∗a2>1,(42)u3∗=0,ifλ3−λ4Iu∗a3<0,u3∗,if0≤λ3−λ4Iu∗a3≤1,1,ifλ3−λ4Iu∗a3>1,u4∗=0,ifλ1S∗a4<0,u4∗,if0≤λ1S∗a4≤1,1,ifλ1S∗a4>1.
More precisely, Equation (42) can be written in the form of Equation (40). The second derivative of the Hamiltonian (H) with respect to u is positive definite. That is, ∂2Hx,u,λ/∂ui2>0. Therefore, the control u∗ is then a minimizer.


## 4.1. Existence of the Optimal Control

Theorem 6.
ConsiderJu subject to system Equation (24) with S0,I0,Iu0,R0≥0,0,0,0, then there exists an optimal control u∗ and corresponding x∗ that minimizes Ju.We need to check the following assumptions based on [31] to prove the existence of an optimal control.
(i)
The set of controls and corresponding state variables are nonempty(ii)
The control setU is convex and closed(iii)
The right hand side of the system Equation (24) is bounded by a linear function in the state and controlIntegrand of the system Equation (26) is convex on U and is bounded below by C1u12+u22+u32+u42μ/2−C2, where C1>0,C2>0 and μ>1.Proof.
(i)
The control variableU is nonempty set of measurable function on the finite interval time 0,T in the real numbers(ii)
Clearly, the control setU is convex and closed by definition(iii)
The right hand side of system Equation (24) is, by definition, continuous on the finite time interval 0,T, and can be written as a linear function of the control U with coefficients depending on time and state. In addition, all the state and control variables are bounded on 0,T
To verify this condition, letλ∈0,1,v=v1,v2,v3,v4∈U, and w=w1,w2,w3,w4∈U.(28)fx,1−λv+λw−1−λfx,v+λfx,w=a121−λ2v12+2λ1−λv1w1+λ2w12+a221−λ2v22+2λ1−λv2w2+λ2w22+a321−λ2v32+2λ1−λv3w3+λ2w32+a421−λ2v42+2λ1−λv4w4+λ2w42−1−λa12v12+a22v22+a32v32+a42v42−λa12w12+a22w22+a32w32+a42w42=λ2−λa12v1−w12+a22v2−w22+a32v3−w32+a42v4−w42=aλ2−λ2v−w2≤0.
It follows that(29)fx,1−λv+λw≤1−λfx,v+λfx,w.(30)fx,u=AI+BIu+a1u122+a2u222+a3u322+a4u422≥C1u12+u22+u32+u42μ/2−C2,where C1=Mina1,a2,a3,a4,μ=2 and C2>0. Thus, the condition is justified.


## 4.2. Characterization of the Optimal Controls

The characterization of the optimal control is based on the Pontryagin’s maximum principle [32]. This principle converts the optimal control problem into the problem of minimizing pointwise a Hamiltonian with respect to u. Then, the Hamiltonian (H) is given, for all t∈0,T, by
(31)Hxt,ut,λt=AI+BIu+a1u122+a2u222+a3u322+a4u422+λ1Λ−βI+u1+u4+μS+λ2βSI−α+u2+μI+λ3αI−δ+u3+μIu+λ4u2I+u3Iu+u1S−μR.Following the approach [30], if the control u∗ and corresponding state x∗ are the optimal control pair, then necessarily there exists a nonzero adjoint vector function λt fulfilling the following equality.
(32)dxtdt=∂Hxt,ut,λt∂λi,controlsystem,(33)dλtdt=−∂Hxt,ut,λt∂xi,adjointcondition,(34)∂Hxt,ut,λt∂ui=0,optimalitycondition,(35)λT=0,transversalitycondition.Finally, it follows that(36)u∗=0,if∂H∂ui<0,0≤ui∗≤1,if∂H∂ui=0,1,if∂H∂ui>0.Theorem 7 (see [29–31]).
Consider the optimal controlu∗t and corresponding trajectory x∗t with t∈0,T for the system Equation (24), then it is necessary that there exist nontrivial adjoint vector function λt satisfying
(37)dλ1dt=βI+u1+u4+μλ1−βIλ2−u1λ4,dλ2dt=βSλ1−βS−α+μ+u2λ2−αλ3−u2λ4−A,dλ3dt=δ+μ+u3λ3−u3λ4−B,dλ4dt=μλ4,λi=0,i=1,2,3,4.
Hence,(38)Hxt,u∗t,λt=minHxt,ut,λt:u∈U.
Furthermore, the compact form of the optimal control,u∗, is given by
(39)u1∗=max0,min1,λ1−λ4S∗a1,u2∗=max0,min1,λ2−λ4I∗a2,(40)u3∗=max0,min1,λ3−λ4Iu∗a3,u4∗=max0,min1,λ1S∗a4.Proof.
We find the adjoint system Equation (37) by differentiating the Hamiltonian (H) in Equation (31) with respect to x=S,I,Iu,R. More precisely, dλt/dt=−∂Hxt,ut,λt/∂xi. We assume that the final state xT is free. Hence, the transversality condition holds that λi=0,i=1,2,3,4. Moreover, the optimal controls, u1∗,u2∗,u3∗,andu4∗, are obtained based on the optimality condition ∂Hx,u,λ/∂ui=0.
Thus, the controls are become as follows:(41)u1∗=0,ifλ1−λ4S∗a1<0,u1∗,if0≤λ1−λ4S∗a1≤1,1,ifλ1−λ4S∗a1>1,u2∗=0,ifλ2−λ4I∗a2<0,u2∗,if0≤λ2−λ4I∗a2≤1,1,ifλ2−λ4I∗a2>1,(42)u3∗=0,ifλ3−λ4Iu∗a3<0,u3∗,if0≤λ3−λ4Iu∗a3≤1,1,ifλ3−λ4Iu∗a3>1,u4∗=0,ifλ1S∗a4<0,u4∗,if0≤λ1S∗a4≤1,1,ifλ1S∗a4>1.
More precisely, Equation (42) can be written in the form of Equation (40). The second derivative of the Hamiltonian (H) with respect to u is positive definite. That is, ∂2Hx,u,λ/∂ui2>0. Therefore, the control u∗ is then a minimizer.


## 5. Numerical Results and Discussions

In this case, we carried out numerical simulation of the model to support our analytical results using reasonable parameter values of the system illustrated in Table2. In Figure 2, numerical simulation of cases without the implementation of any optimal control strategy against time can be seen. In this regard, the number of cases increases initially and then slowly drops to its endemic equilibrium. From Figure 2, it is evident that the infectious individuals, I, reach its maximum value 2127 at the approximate period of time t=5.28 days, while the infectious population under intensive care unit, Iu, attains its maximum value 1641 at an approximate time t=15 days from 2200 population.Table 2
The initial values for the state variables and parameter values of system model.
ParameterValueSourceState variableInitial valueSourceΛ100AssumedS03000Assumedβ0.000218AssumedI0600Assumedν0.0001AssumedIu0400Assumedμ4.7788×10–5[33]R0300Assumedα0.1Assumedγ10.03571[34]γ20.04255[35]δ0.0212AssumedFigure 2
Number of cases without any optimal control strategy versus time.As can be observed in Figures3(a) and 3(b), numerical simulation of the infected population has been carried out with different initial values for β=0.0000218&ν=0.1 (hence, R0=0.1530<1). The other parameters are as indicated in Table 2. These numerical results illustrate the global stability of the disease-free equilibrium of the system Equation (3). Moreover, Figures 3(c) and 3(d) present the time series of system Equation (3) for parameters illustrated in Table 2 and, hence, R0=1.5296>1. From this numerical simulation, one can observe that the solution of system Equation (3) ultimately converges to its endemic equilibrium. This is an evident that our numerical results are in a good agreement with the analytical results.Figure 3
Model simulation with different initial conditions.
(a)
The infectedI against time for R0<1(b)
The infectedIu against time for R0<1(c)
The infectedI against time for R0>1(d)
The infectedIu against time for R0>1


### 5.1. The Effect of Optimal Control Strategies

In this section, we examine the results of numerical simulations with the implementation of optimal control strategies for systems Equation (24).We employ the forward-backward sweep method to solve the optimality system. The procedure can be summarized as follows: first, we solve the control system Equation (32) with an initial guess value using the forward fourth-order Runge-Kutta method. By employing the backward fourth-order Runge-Kutta method, we therefore solved the adjoint system Equation (37) using the current iteration solution of the state variables. The controls are updated by using a convex combination of the previous controls and the values from Equation (42). The process continues till the results at the consecutive iterations are too enough.The following are chosen as parameter values for the optimal control simulation:Λ=100,β=0.000218,μ=4.7788×10−5,α=0.01,δ=0.0212. Moreover, the assumed weight factors are A=10,B=5,a1=5,a2=10,a3=10,a4=5. Among the I (infected) and Iu (ICU), the number of infected humans is considered the most important ones because it is assumed in the model formulation that disease transmission is due to the contact between the susceptible and the infected groups. Therefore, the highest attention A=10 is assumed to be to I and relatively less (B=5) to Iu. The cost of controlling the disease is relatively smaller but still important. For example, health education through creating awareness and vaccination is comparatively cheaper and assumed to be a1=a4=5, whereas treatment is relatively more expensive; hence, set a2=a3=10. These values may differ depending on different regions and nations. On the assumption that the maximum rate of intervention per day is 100%, we set u1max=u2max=u3max=u4max=1. All the optimal control simulations are carried out by using the initial state variables given in Table 2.Accordingly, we study numerically to show the impact of different control strategies on the control of the spread of the disease.(i)
Strategy A. Implementing all controls u1,u2,u3, and u4(ii)
Strategy B. Implementing controls u1,u3, and u4(iii)
Strategy C. Implementing controls u1, and u3(iv)
Strategy D. Implementing controls u1 and u4(v)
Strategy A. Implementing controls u2 and u3


#### 5.1.1. Strategy A: Implementing All Controlsu1, u2, u3, and u4

In this case, we use all the control strategies to optimize the objective function in Equation (26). There is a significant difference in number of cases in simulation of the system with and without controls. It can be observed from Figure 2 that the number of cases is plotted without any controls. In this case, the number of infections showed an increasing alarm and slowly adjusted to their endemic equilibrium. On the other hand, Figures 4(a) and 4(b) indicate that the number of cases I and Iu quickly decreased and lead to zero, respectively, at times t=6.5 and t=8 days. Specifically, the total number of infectious, I and Iu, at the end of 120 days become 743 and 1565, respectively, without controls, and zero with implementing all control strategies. This result shows the significant effect of the optimal control measures on the dynamics of the disease.Figure 4
Number of cases,I and Iu, with all control strategies.
(a)
I with and without controls(b)
Iu with and without controlsFrom the control profile Figure5(a), the control effort u3 maintained its maximum level till time t=16.7 days and gradually dropped to its lower bound in the remaining time. To minimize the burden of the disease, the optimal control measures u1 and u2 should be kept at their maximum level for about the first 18.5 and 19 days, respectively, while the control effort u4 remains at its maximum rate for the whole period.Figure 5
Plot of optimal control profiles.
(a)
Optimal control profiles whenui≠0,i=1,2,3,4(b)
Optimal control profiles whenui≠0,i=1,3,4(c)
Optimal control profiles whenui≠0,i=1,3(d)
Optimal control profiles whenui≠0,i=1,4(e)
Optimal control profiles whenui≠0,i=2,3


#### 5.1.2. Strategy B: Implementing Controlsu1, u3, and u4

This is the case where the vaccination controlu1, the treatment control u3, and the public health education control u4 are applied to minimize the infected population while we set the treatment control u2 to zero. One can observe from Figure 6(a) that due to the implementation of control measures, an increased number is observed for the first one day in the infectious individuals I and then slowly decreased to 135 for the remaining periods, while the infected population under intensive care unit Iu gradually declined and reaches to its lower value 16 at time t=20 days Figure 6(b). In Figures 6(a) and 6(b), it can be observed that the number of cases is reduced when optimal controls are implemented compared to uncontrolled ones. The optimal controls for vaccination u1, treatment u3, and public health education u4 are illustrated in Figure 5(b).Figure 6
Number of cases,I and Iu, with optimal control profiles u1,u3, and u4.
(a)
I with and without interventions(b)
Iu with and without controls


#### 5.1.3. Strategy C: Implementing Controlsu1 and u3

To optimize our objective function, we used the optimal control vaccinationu1 and the treatment control u3, while the treatment control u2 and the public health education control u4 are set to zero. It can be observed from Figures 7(a) and 7(b) that there is a significant difference in the number of cases with and without controls. That means, the number of infectious individuals I with controls shows an initial increment and gradually drops to 208 at the time t=20 days, while the infected population under the intensive care unit Iu revealed quickly declined at the initial time and tends to 23 at the end of the period, whereas the number of cases I and Iu without controls shows a rapid increased at the initial time and slowly decreased until they adjusted to their stable endemic equilibrium. More importantly, at the end of 120 days, the total number of cases I+Iu is 2308 without controls and 231 with implementing the controls.Figure 7
The number of cases,I and Iu, with control strategies u1 and u3.
(a)
The number of cases,I, with and without controls(b)
The number of cases,Iu, with and without controlsTo reduce the burden of infections in the community, we implemented the controlsu1 and u3 illustrated in Figure 5(c). The results in Figure 5(c) show that the vaccination control u1 stayed in its maximum level for the first 19 days and afterwards dropped to its lower bound for the remaining periods, while the treatment rate u3 is at its maximum rate until the period of 19.5 days and be declined to its lower bound for about the period of half a day.


#### 5.1.4. Strategy D: Implementing Controlsu1 and u4

In this case, we implemented the vaccination controlu1 and the public health education control u4 in order to minimize the number of infected people in the absence of treatment controls u2 and u3. The corresponding simulation results are presented in Figures 8(a) and 8(b). From this, we notice that the infected people I increases and reaches its maximum peak 716 at time of one day, while the infectious individuals under the intensive care unit Iu attain its maximum peak 832 during the whole periods. However, the maximum peak values of the cases I and Iu without implemented controls become 2127 at an approximate time t=5.4 days and 1641 at an approximate time t=15 days, respectively. This shows that there is a significant difference in the number of cases with and without controls.Figure 8
Changes in number of cases,I and Iu, with the optimal control profiles u1 and u4.
(a)
Changes in number of cases,I, with and without controls(b)
Changes in number of cases,Iu, with and without controlsFigure5(d) shows the control profiles of the vaccination u1 and the public health education u4 in which the optimal controls u1 and u4 are maintained at their upper bound for about the first 19.4 and the whole periods, respectively. Moreover, the control u1 drops rapidly to its lower bound for the remaining time.


#### 5.1.5. Strategy E: Implementing Controlsu2 and u3

We use the two treatment controls,u2 and u3 with u1=0 and u4=0, to optimize the objective functional J found in Equation (26). We observed from Figure 9(a) that the number of infectious individuals I shows a sharp declined and then tends to zero near to the time t=16 days, while the number of infectious population with the intensive care unit Iu leads to zero in the approximate time t=20 days as can be seen in Figure 9(b).Figure 9
Numerical simulation for the number of cases,I and Iu, with the optimal control profiles u2 and u3.
(a)
Time series for the number of cases,I, with and without controls(b)
Time series for the number of cases,Iu, with and without controlsFigure5(a) illustrates the controls u2 and u3 against time. To minimize the number of cases, the optimal control u2 is kept at its upper bound for about 18 days and then gradually dropped to its lower bound, whereas the optimal control u3 is kept at its upper bound for about 19.8 days and then sharply dropped to its lower bound. The treatment controls u2 and u3 in this strategy have to be prolonged as compared to in strategy A before dropping to their lower bounds (see Figures 5(a) and 5(e)).In conclusion, we observed from the optimal control strategies that the combined all controls and the combined treatment controls (u2 and u3) show the most significant effect in eradicating the burden of COVID-19 in society relative to the other combinations. The possible implication is that the disease can be eradicated by implementing the combination of all controls and or the combination of the treatments u2 and u3. However, implementation of the intervention measures u1 and u4 or their combinations can provide an important contribution in reducing the peak values in a number of cases, as can be seen in Figures 8(a) and 8(b).


## 5.1. The Effect of Optimal Control Strategies

In this section, we examine the results of numerical simulations with the implementation of optimal control strategies for systems Equation (24).We employ the forward-backward sweep method to solve the optimality system. The procedure can be summarized as follows: first, we solve the control system Equation (32) with an initial guess value using the forward fourth-order Runge-Kutta method. By employing the backward fourth-order Runge-Kutta method, we therefore solved the adjoint system Equation (37) using the current iteration solution of the state variables. The controls are updated by using a convex combination of the previous controls and the values from Equation (42). The process continues till the results at the consecutive iterations are too enough.The following are chosen as parameter values for the optimal control simulation:Λ=100,β=0.000218,μ=4.7788×10−5,α=0.01,δ=0.0212. Moreover, the assumed weight factors are A=10,B=5,a1=5,a2=10,a3=10,a4=5. Among the I (infected) and Iu (ICU), the number of infected humans is considered the most important ones because it is assumed in the model formulation that disease transmission is due to the contact between the susceptible and the infected groups. Therefore, the highest attention A=10 is assumed to be to I and relatively less (B=5) to Iu. The cost of controlling the disease is relatively smaller but still important. For example, health education through creating awareness and vaccination is comparatively cheaper and assumed to be a1=a4=5, whereas treatment is relatively more expensive; hence, set a2=a3=10. These values may differ depending on different regions and nations. On the assumption that the maximum rate of intervention per day is 100%, we set u1max=u2max=u3max=u4max=1. All the optimal control simulations are carried out by using the initial state variables given in Table 2.Accordingly, we study numerically to show the impact of different control strategies on the control of the spread of the disease.(i)
Strategy A. Implementing all controls u1,u2,u3, and u4(ii)
Strategy B. Implementing controls u1,u3, and u4(iii)
Strategy C. Implementing controls u1, and u3(iv)
Strategy D. Implementing controls u1 and u4(v)
Strategy A. Implementing controls u2 and u3


### 5.1.1. Strategy A: Implementing All Controlsu1, u2, u3, and u4

In this case, we use all the control strategies to optimize the objective function in Equation (26). There is a significant difference in number of cases in simulation of the system with and without controls. It can be observed from Figure 2 that the number of cases is plotted without any controls. In this case, the number of infections showed an increasing alarm and slowly adjusted to their endemic equilibrium. On the other hand, Figures 4(a) and 4(b) indicate that the number of cases I and Iu quickly decreased and lead to zero, respectively, at times t=6.5 and t=8 days. Specifically, the total number of infectious, I and Iu, at the end of 120 days become 743 and 1565, respectively, without controls, and zero with implementing all control strategies. This result shows the significant effect of the optimal control measures on the dynamics of the disease.Figure 4
Number of cases,I and Iu, with all control strategies.
(a)
I with and without controls(b)
Iu with and without controlsFrom the control profile Figure5(a), the control effort u3 maintained its maximum level till time t=16.7 days and gradually dropped to its lower bound in the remaining time. To minimize the burden of the disease, the optimal control measures u1 and u2 should be kept at their maximum level for about the first 18.5 and 19 days, respectively, while the control effort u4 remains at its maximum rate for the whole period.Figure 5
Plot of optimal control profiles.
(a)
Optimal control profiles whenui≠0,i=1,2,3,4(b)
Optimal control profiles whenui≠0,i=1,3,4(c)
Optimal control profiles whenui≠0,i=1,3(d)
Optimal control profiles whenui≠0,i=1,4(e)
Optimal control profiles whenui≠0,i=2,3


### 5.1.2. Strategy B: Implementing Controlsu1, u3, and u4

This is the case where the vaccination controlu1, the treatment control u3, and the public health education control u4 are applied to minimize the infected population while we set the treatment control u2 to zero. One can observe from Figure 6(a) that due to the implementation of control measures, an increased number is observed for the first one day in the infectious individuals I and then slowly decreased to 135 for the remaining periods, while the infected population under intensive care unit Iu gradually declined and reaches to its lower value 16 at time t=20 days Figure 6(b). In Figures 6(a) and 6(b), it can be observed that the number of cases is reduced when optimal controls are implemented compared to uncontrolled ones. The optimal controls for vaccination u1, treatment u3, and public health education u4 are illustrated in Figure 5(b).Figure 6
Number of cases,I and Iu, with optimal control profiles u1,u3, and u4.
(a)
I with and without interventions(b)
Iu with and without controls


### 5.1.3. Strategy C: Implementing Controlsu1 and u3

To optimize our objective function, we used the optimal control vaccinationu1 and the treatment control u3, while the treatment control u2 and the public health education control u4 are set to zero. It can be observed from Figures 7(a) and 7(b) that there is a significant difference in the number of cases with and without controls. That means, the number of infectious individuals I with controls shows an initial increment and gradually drops to 208 at the time t=20 days, while the infected population under the intensive care unit Iu revealed quickly declined at the initial time and tends to 23 at the end of the period, whereas the number of cases I and Iu without controls shows a rapid increased at the initial time and slowly decreased until they adjusted to their stable endemic equilibrium. More importantly, at the end of 120 days, the total number of cases I+Iu is 2308 without controls and 231 with implementing the controls.Figure 7
The number of cases,I and Iu, with control strategies u1 and u3.
(a)
The number of cases,I, with and without controls(b)
The number of cases,Iu, with and without controlsTo reduce the burden of infections in the community, we implemented the controlsu1 and u3 illustrated in Figure 5(c). The results in Figure 5(c) show that the vaccination control u1 stayed in its maximum level for the first 19 days and afterwards dropped to its lower bound for the remaining periods, while the treatment rate u3 is at its maximum rate until the period of 19.5 days and be declined to its lower bound for about the period of half a day.


### 5.1.4. Strategy D: Implementing Controlsu1 and u4

In this case, we implemented the vaccination controlu1 and the public health education control u4 in order to minimize the number of infected people in the absence of treatment controls u2 and u3. The corresponding simulation results are presented in Figures 8(a) and 8(b). From this, we notice that the infected people I increases and reaches its maximum peak 716 at time of one day, while the infectious individuals under the intensive care unit Iu attain its maximum peak 832 during the whole periods. However, the maximum peak values of the cases I and Iu without implemented controls become 2127 at an approximate time t=5.4 days and 1641 at an approximate time t=15 days, respectively. This shows that there is a significant difference in the number of cases with and without controls.Figure 8
Changes in number of cases,I and Iu, with the optimal control profiles u1 and u4.
(a)
Changes in number of cases,I, with and without controls(b)
Changes in number of cases,Iu, with and without controlsFigure5(d) shows the control profiles of the vaccination u1 and the public health education u4 in which the optimal controls u1 and u4 are maintained at their upper bound for about the first 19.4 and the whole periods, respectively. Moreover, the control u1 drops rapidly to its lower bound for the remaining time.


### 5.1.5. Strategy E: Implementing Controlsu2 and u3

We use the two treatment controls,u2 and u3 with u1=0 and u4=0, to optimize the objective functional J found in Equation (26). We observed from Figure 9(a) that the number of infectious individuals I shows a sharp declined and then tends to zero near to the time t=16 days, while the number of infectious population with the intensive care unit Iu leads to zero in the approximate time t=20 days as can be seen in Figure 9(b).Figure 9
Numerical simulation for the number of cases,I and Iu, with the optimal control profiles u2 and u3.
(a)
Time series for the number of cases,I, with and without controls(b)
Time series for the number of cases,Iu, with and without controlsFigure5(a) illustrates the controls u2 and u3 against time. To minimize the number of cases, the optimal control u2 is kept at its upper bound for about 18 days and then gradually dropped to its lower bound, whereas the optimal control u3 is kept at its upper bound for about 19.8 days and then sharply dropped to its lower bound. The treatment controls u2 and u3 in this strategy have to be prolonged as compared to in strategy A before dropping to their lower bounds (see Figures 5(a) and 5(e)).In conclusion, we observed from the optimal control strategies that the combined all controls and the combined treatment controls (u2 and u3) show the most significant effect in eradicating the burden of COVID-19 in society relative to the other combinations. The possible implication is that the disease can be eradicated by implementing the combination of all controls and or the combination of the treatments u2 and u3. However, implementation of the intervention measures u1 and u4 or their combinations can provide an important contribution in reducing the peak values in a number of cases, as can be seen in Figures 8(a) and 8(b).


## 5.1.1. Strategy A: Implementing All Controlsu1, u2, u3, and u4

In this case, we use all the control strategies to optimize the objective function in Equation (26). There is a significant difference in number of cases in simulation of the system with and without controls. It can be observed from Figure 2 that the number of cases is plotted without any controls. In this case, the number of infections showed an increasing alarm and slowly adjusted to their endemic equilibrium. On the other hand, Figures 4(a) and 4(b) indicate that the number of cases I and Iu quickly decreased and lead to zero, respectively, at times t=6.5 and t=8 days. Specifically, the total number of infectious, I and Iu, at the end of 120 days become 743 and 1565, respectively, without controls, and zero with implementing all control strategies. This result shows the significant effect of the optimal control measures on the dynamics of the disease.Figure 4
Number of cases,I and Iu, with all control strategies.
(a)
I with and without controls(b)
Iu with and without controlsFrom the control profile Figure5(a), the control effort u3 maintained its maximum level till time t=16.7 days and gradually dropped to its lower bound in the remaining time. To minimize the burden of the disease, the optimal control measures u1 and u2 should be kept at their maximum level for about the first 18.5 and 19 days, respectively, while the control effort u4 remains at its maximum rate for the whole period.Figure 5
Plot of optimal control profiles.
(a)
Optimal control profiles whenui≠0,i=1,2,3,4(b)
Optimal control profiles whenui≠0,i=1,3,4(c)
Optimal control profiles whenui≠0,i=1,3(d)
Optimal control profiles whenui≠0,i=1,4(e)
Optimal control profiles whenui≠0,i=2,3


## 5.1.2. Strategy B: Implementing Controlsu1, u3, and u4

This is the case where the vaccination controlu1, the treatment control u3, and the public health education control u4 are applied to minimize the infected population while we set the treatment control u2 to zero. One can observe from Figure 6(a) that due to the implementation of control measures, an increased number is observed for the first one day in the infectious individuals I and then slowly decreased to 135 for the remaining periods, while the infected population under intensive care unit Iu gradually declined and reaches to its lower value 16 at time t=20 days Figure 6(b). In Figures 6(a) and 6(b), it can be observed that the number of cases is reduced when optimal controls are implemented compared to uncontrolled ones. The optimal controls for vaccination u1, treatment u3, and public health education u4 are illustrated in Figure 5(b).Figure 6
Number of cases,I and Iu, with optimal control profiles u1,u3, and u4.
(a)
I with and without interventions(b)
Iu with and without controls


## 5.1.3. Strategy C: Implementing Controlsu1 and u3

To optimize our objective function, we used the optimal control vaccinationu1 and the treatment control u3, while the treatment control u2 and the public health education control u4 are set to zero. It can be observed from Figures 7(a) and 7(b) that there is a significant difference in the number of cases with and without controls. That means, the number of infectious individuals I with controls shows an initial increment and gradually drops to 208 at the time t=20 days, while the infected population under the intensive care unit Iu revealed quickly declined at the initial time and tends to 23 at the end of the period, whereas the number of cases I and Iu without controls shows a rapid increased at the initial time and slowly decreased until they adjusted to their stable endemic equilibrium. More importantly, at the end of 120 days, the total number of cases I+Iu is 2308 without controls and 231 with implementing the controls.Figure 7
The number of cases,I and Iu, with control strategies u1 and u3.
(a)
The number of cases,I, with and without controls(b)
The number of cases,Iu, with and without controlsTo reduce the burden of infections in the community, we implemented the controlsu1 and u3 illustrated in Figure 5(c). The results in Figure 5(c) show that the vaccination control u1 stayed in its maximum level for the first 19 days and afterwards dropped to its lower bound for the remaining periods, while the treatment rate u3 is at its maximum rate until the period of 19.5 days and be declined to its lower bound for about the period of half a day.


## 5.1.4. Strategy D: Implementing Controlsu1 and u4

In this case, we implemented the vaccination controlu1 and the public health education control u4 in order to minimize the number of infected people in the absence of treatment controls u2 and u3. The corresponding simulation results are presented in Figures 8(a) and 8(b). From this, we notice that the infected people I increases and reaches its maximum peak 716 at time of one day, while the infectious individuals under the intensive care unit Iu attain its maximum peak 832 during the whole periods. However, the maximum peak values of the cases I and Iu without implemented controls become 2127 at an approximate time t=5.4 days and 1641 at an approximate time t=15 days, respectively. This shows that there is a significant difference in the number of cases with and without controls.Figure 8
Changes in number of cases,I and Iu, with the optimal control profiles u1 and u4.
(a)
Changes in number of cases,I, with and without controls(b)
Changes in number of cases,Iu, with and without controlsFigure5(d) shows the control profiles of the vaccination u1 and the public health education u4 in which the optimal controls u1 and u4 are maintained at their upper bound for about the first 19.4 and the whole periods, respectively. Moreover, the control u1 drops rapidly to its lower bound for the remaining time.


## 5.1.5. Strategy E: Implementing Controlsu2 and u3

We use the two treatment controls,u2 and u3 with u1=0 and u4=0, to optimize the objective functional J found in Equation (26). We observed from Figure 9(a) that the number of infectious individuals I shows a sharp declined and then tends to zero near to the time t=16 days, while the number of infectious population with the intensive care unit Iu leads to zero in the approximate time t=20 days as can be seen in Figure 9(b).Figure 9
Numerical simulation for the number of cases,I and Iu, with the optimal control profiles u2 and u3.
(a)
Time series for the number of cases,I, with and without controls(b)
Time series for the number of cases,Iu, with and without controlsFigure5(a) illustrates the controls u2 and u3 against time. To minimize the number of cases, the optimal control u2 is kept at its upper bound for about 18 days and then gradually dropped to its lower bound, whereas the optimal control u3 is kept at its upper bound for about 19.8 days and then sharply dropped to its lower bound. The treatment controls u2 and u3 in this strategy have to be prolonged as compared to in strategy A before dropping to their lower bounds (see Figures 5(a) and 5(e)).In conclusion, we observed from the optimal control strategies that the combined all controls and the combined treatment controls (u2 and u3) show the most significant effect in eradicating the burden of COVID-19 in society relative to the other combinations. The possible implication is that the disease can be eradicated by implementing the combination of all controls and or the combination of the treatments u2 and u3. However, implementation of the intervention measures u1 and u4 or their combinations can provide an important contribution in reducing the peak values in a number of cases, as can be seen in Figures 8(a) and 8(b).


## 6. Conclusion

In this study, we formulated and analyzed a mathematical model with vaccination (v) as a constant control that describes the dynamics of COVID-19 disease. We analyzed the model for the existence and global stability of disease-free and endemic equilibria. For global stability, we used Lyapunov function methods. In particular, the global stability analysis of the endemic equilibrium is done by constructing a logarithmic Lyapunov function in the SI plane, which is similar to the works of Beretta et al., Korobeinikov et al., and Vargas-De Leon [36–38]. The system formulated is further extended into an optimal control problem by introducing continuous controls: vaccination, treatments, and public health education into the system to minimize the burden of the disease in the community and the associated costs. We verified analytically the existence of the optimal controls and characterized them by employing Pontryagin’s maximum principle. We also numerically conducted comparative studies by designing different combinations of control strategies to propose better intervention measures in the reduction and eradication of the disease.The following can be highlighted from the analysis in this study:(i)
It is found from the mathematical analysis that both the disease-free and endemic equilibria are globally asymptotically stable by constructing suitable Lyapunov functions(ii)
We found also from the local sensitivity analysis that the most sensitive parameters are the recruitment rateΛ, the transmission rate β, and the vaccination rate ν. It is moreover found out that the reproduction number R0 is independent of the parameter recovery rate of the ICU γ2 and the death rate δ of the ICU. Thus, the spread of the disease does not depend on the increase or decrease of these parametersThis work has limitations. The model formulated is deterministic and based on homogeneous mixing with the assumption that all susceptible individuals have the same chance of being infected with the disease if they come into contact with the infected ones [2, 5, 7, 30, 39]. Another limitation in this study is the assumption that the ICU is not the source of infection. But, it is vital to consider an epidemic model with the ICU as the source of infection in addition to the outpatients. It is also important that the model proposed be validated with real data. It is assumed in all the above cases that the disease emerges in a large population. However, if the population is small, susceptible individuals can possibly contact the infectious at random, and reporting the number of contacts as a function of the population density is difficult [40]. In this case, a stochastic model is preferred. A mathematical model for COVID-19 disease that fills these gaps will be considered part of our future work.



---

*Source: 8636530-2022-11-09.xml*