Tinaw Tilahun Asmamaw1, Kiros Gebrearegawi Kebedow2,∗ 1Department of Mathematics, Mezan Tepi University, Ethiopia 2Department of Mathematics, Hawassa University, Ethiopia KEYWORDS COVID-19; Protection; Stability analysis; Forward bifurcation; Sensitivity analysis; Optimal control ABSTRACT ∗ East Afr. J. Biophys. Comput. Sci. (2024), Issue 5, No. 2, 58-86 58 In this paper, we propose a mathematical model to investigate coronavirus diseases (COVID-19) transmission in the presence of protected and hos- pitalized classes. We establish that the solution of the dynamical system remains positive and bounded. We compute the disease free equilibrium point and analyze the stability behavior of the steady state solutions.We determine the basic reproduction number (R0)and demonstrate that the disease fades away when R0 < 1but persists in the population when R0 > 1. The center manifold theory is used to assess the local stability of the endemic equilibrium. The model demonstrates a forward bifurcation, and a sensitivity analysis is conducted. The sensitivity analysis reveals that R0 is highly influenced by the protection rate, highlighting the necessity of maintaining a high level of protection along with hospitalization to effectively control the disease. We develop optimal strategies for protection and hospitalization. The characterization of the optimal control is derived using Pontryagin’s Maximum Principle. Numerical results for the dynamics of the COVID-19 outbreak and its optimal control show that a combination of protection and hospitalization is the most effective strategy for reducing the spread of COVID-19 within the population.. portation Lemecha Obsu and Feyissa Balcha (2020). The novel coronavirus- now referred to as COVID-19 is caused by severe acute respira- tory syndrome (SARS-CoV-2) and consists of single-stranded ribonucleic acid (RNA) struc- ture Sohrabi et al. (2020). The novel coronavirus is mainly spread 1. INTRODUCTION Coronavirus disease 2019 (COVID-19) is an infectious disease caused by a newly discov- ered coronavirus Gurmu et al. (2020). The new virus was first appeared late December 2019 in the Chinese city of Wuhan and eventually in- vaded the world due to fast modern air trans- A Mathematical Model for the Transmission Dynamics of COVID-19 Pandemic Considering Protected and Hospitalized with Optimal Control Corresponding author Email address: kirosg@hu.edu.et +251926528484 https://dx.doi.org/10.4314/eajbcs.v5i2.5S Research article East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 59 from person to person, through respiratory droplets, the spread is more likely when peo- ple are within 6 feet of each other Toquero (2020). There is no known curing medicine to combat the COVID-19 pandemic. Standard recommendations by the WHO to prevent the spread of COVID-19 include frequent cleaning of hands using soap or alcohol-based sanitizer, covering the nose and mouth with a flexed el- bow or disposable tissue when coughing and sneezing and avoiding close contact with any- one that has a fever and cough Toquero (2020). Mathematical modeling in epidemiology helps to understand the fundamental mechanisms that drive the spread of disease, while also offering insights into potential control strategies.The model formulation process clarifies the assumptions, variables, and parameters involved. Furthermore, models offer conceptual insights such as thresholds, basic reproduction numbers, contact numbers, and replacement numbers.Mathematical models and computer simulations are useful experimental tools for building and testing theories, assessing quan-titative conjectures, answering specific ques-tions, determining sensitivities to changes in parameter values, and estimating parameters from data Brauer et al. (2019). The first m athematical m odel w as devel- oped by Daniel Bernoulli Allen et al. (2008) on pandemic of smallpox by introducing two systems of ordinary differential e quations. He assumed that recovery from infection confers immunity (no re-infection). He also assumed that the probability that infected individuals for the first time die does not depend on those who survive from the infection. He showed that inoculation was advantageous if the as- sociated risk of dying was less than 11%. A number of compartmental models have been proposed and analyzed for the COVID- 19 outbreak in different c ountries. I n par- ticular, Yang and Wang Yang and Wang (2020) proposed a mathematical model for COVID-19 incorporating multiple transmis- sion pathways, including both human-to- human and environment-to-human transmis- sion routes. Global stability was analysed us- ing Lyapunov function. The authors employed a bilinear incidence rate based on the law of mass action and fitted the model with the data of Wuhan city of China and estimated the re- production number. In 2020, Haileyesus and Getachew Alem- neh and Telahun (2020) proposed a conceptual SEIR model to study the pandemic COVID-19 transmission in Ethiopia. Global stability was analysed using Lyapunov function. Additionally, they incorporated time- dependent controls into the basic model and extended it to an optimal control framework for the disease. An optimal control problem was formulated and analyzed using Pontryagin’s Maximum Principle. However, none of the authors cited here are considering protected and hospital-ized classes with optimal controls. Thus we are improved the existed model, by including the protected and hospitalized classes with op- timal controls. The remainder of the paper is organized as follows. In Section 2 we provided the de- scription of the problems and its mathemati- cal model formulation. Section 3 gives details the qualitative analysis of the model. The nu- merical analysis of our model is presented in Section 4, furthermore we present the exten- sions of the model to optimal controls and its numerical simulations 5 and 6. Finally, Section 7 concludes the paper. 2. MATHEMATICAL MODEL FOR- MULATION 2.1. The modified mathematical model In this subsection, we modify the exist- ing SEIR model to study the transmission dy- namics of COVID-19 infection in a population. The model is a modification of what is pre- sented in Alemneh and Telahun (2020). In the present model we extended SEIR model by in- cluding protected and hospitalized classes. 2.2. Model formulation We formulate a mathematical model to investigate coronavirus diseases (COVID-19) transmission in the presence of protected and hospitalized classes. The model divides the to- tal population into six sub-classes according to their disease status. Susceptible (S), protected (P), exposed (E), infected (I), hospitalized (H) and recovered (R). The following assumptions have been used in the formulation of the model: 1. The population under study is heteroge- neous and varying with time. 2. All recruited human population is sus- ceptible. 3. Susceptible individuals who keeping so- cial distancing, using an alcohol-based hand sanitizer and wearing face masks progress into protected class. 4. Protected individuals cannot acquire in- fection of COVID-19 disease due to proper use of an alcohol-based hand san- itizer, wearing face masks and keeping social distancing. 5. The latently infected individuals (ex- posed) could infect other people with a higher probability than people in the in- fected class. Since exposed individuals show no symptoms and can easily spread the infection to other people with close contact, often in an unconscious manner. 6. Transmission through human-to-human route is alone considered. Other means of transmission are ignored. 7. We assume that individuals have no per- manent immunity after recovery from the disease, that is the recovered indi- viduals have a chance to be susceptible again. 8 The natural mortality rates are assumed to be the same for all the compartments. 9. All parameters in the model being non- negative. We consider the force of infection λ which is given by λ = β(σ1I + σ2E + σ3H), (1) where β is the effective contact rate, while σ1, σ2 and σ3 are the relative infectiousness pa- rameters associated with the infected, exposed and hospitalized classes respectively. Variables Description of the state variables S Susceptible individuals P Protected individuals, who keeping social distancing, using an alcohol-based hand sanitizer and wearing face masks properly to protect themselves from the virus E Latently infected individuals, who have no symptoms of COVID-19 virus disease but are capable of infecting others I Infected individuals, who have active COVID-19 virus disease and can infect other people H Hospitalized individuals, who are admitted to health care facility or isolated in their home due to virus infection active cases R Recovered individuals The flow chart of the modified model is illustrated in Figure (1). East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 60 Table 1: The state variables and their descriptions. Table 2: Parameters of the modified model and their descriptions. Parameter Description of the parameter Π Recruitment rate of susceptible θ Protection rate of susceptible individuals ν Waning rate of protected individuals to susceptible class β Effective contact rate σ1 Modification parameter for relative infectiousness of infected individuals σ2 Modification parameter for relative infectiousness of exposed individuals σ3 Modification parameter for relative infectiousness of hospitalized individuals δ The exposed progression rate τ Proportion of exposed individuals who join infected class 1− τ The progression from exposed individuals to recovered class α Hospitalization rate of infected individuals µ Natural death rate ρ Disease-induced death rate of infected individuals ξ Disease-induced death rate of hospitalized individuals ε Rate of recovery of the individuals from infected class γ Recovery rate of hospitalized patients ω Waning immunity rate η The proportion of recovered individuals that become susceptible 1− η The progression from recovered individuals to protected class Figure 1: The flow chart for the modified model of COVID-19 pandemic. East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 61 Based on our assumptions and the flow chart (1), the modified model for the transmission dynamics of COVID-19 is given by the following deterministic system of non-linear differential equations: dS dt = Π+ ηωR + νP − β(σ1I + σ2E + σ3H)S − (θ + µ)S, (2a) dP dt = θS + (1− η)ωR− (ν + µ)P, (2b) dE dt = β(σ1I + σ2E + σ3H)S − (δ + µ)E, (2c) dI dt = τδE − (α+ ε+ µ+ ρ)I, (2d) dH dt = αI − (γ + µ+ ξ)H, (2e) dR dt = (1− τ)δE + εI + γH − (ω + µ)R, (2f) with non-negative initial conditions S(0) = S0 > 0, P (0) = P0 ≥ 0, E(0) = E0 ≥ 0, I(0) = I0 ≥ 0, H(0) = H0 ≥ 0 and R(0) = R0 ≥ 0. 3. QUALITATIVE ANALYSIS OF THE MODIFIED MODEL In this section, we present some basic quali- tative properties of the modified model. These analysis include finding t he s et i nside which the model can be sufficiently studied (i .e., the invariant region); local and global stability of equilibrium points of the model (2). 3.1. Well-posedness Since all the functions on the right hand side of the system (2) are continuously differ- entiable. Thus, the existence and uniqueness of the solutions is established by the Picard’s theorem Valcher (2002). Now, we show the positivity and boundedness of solutions. Theorem 3.1. (Positivity) If S(0) > 0, P (0) ≥ 0, E(0) ≥ 0, I(0) ≥ 0, H(0) ≥ 0 and R(0) ≥ 0, then the so- lution (S(t), P (t), E(t), I(t), H(t), R(t)) of the dynamical system (2) is non-negative for all time t ≥ 0. Proof. To show the positivity of the solution of the dynamical system (2), we will perform the proof by using contradiction. We assume that S(t) ≤ 0 for some t ≥ 0, that is there exists small t0 > 0 such that S(t0) = 0, S ′(t0) ≤ 0 and S(t) > 0 for t ∈ [0, t0). Then P (t) ≥ 0, E(t) ≥ 0 and I(t) ≥ 0 for t ∈ [0, t0]. If this be not the case, there exists Option-I: t1 ∈ [0, t0] such that P (t1) = 0, P ′(t1) < 0 and P (t) > 0 for t ∈ [0, t1). Then E(t) ≥ 0 and I(t) ≥ 0 for t ∈ [0, t1]. Option-II: t2 ∈ [0, t1] such that E(t2) = 0, E ′(t2) < 0 and E(t) > 0 for t ∈ [0, t2). Then P (t) ≥ 0 and I(t) ≥ 0 for t ∈ [0, t2]. Option-III: t3 ∈ [0, t2] such that I(t3) = 0, I ′(t3) < 0 and I(t) > 0 for t ∈ [0, t3). Then P (t) ≥ 0 and E(t) ≥ 0 for t ∈ [0, t3]. It follows from equation (2d) that we have I ′(t3) = τδE(t3)− (α+ ε+ µ+ ρ)I(t3). This implies that I ′(t3) = τδE(t3) ≥ 0. This is a contradiction. Integration of equation (2e) leads to H(t) = e−(γ+µ+ξ)t ( H(0) + α ∫ t 0 I(s)e(γ+µ+ξ)sds ) ≥ 0, for t ∈ [0, t2]. East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 62 Then E ′(t2) = β(σ1I(t2)+σ3H(t2))S(t2) ≥ 0. This is a contradiction. Hence H(t) ≥ 0 for every t ∈ [0, t1]. Integration of equation (2f) leads to R(t) = e−(ω+µ)t ( R(0) + ∫ t 0 ((1− τ)δE(s) + εI(s) + γH(s)) e(ω+µ)sds ) ≥ 0, for t ∈ [0, t1]. Then P ′(t1) = θS(t1) + (1 − η)ωR(t1) ≥ 0. This is a contradiction. Hence R(t) ≥ 0 for every t ∈ [0, t0]. Thus S ′(t0) = Π + ηωR(t0) + νP (t0) > 0, but this leads to a contradiction to the assumption that S ′(t0) ≤ 0. Therefore, the solutions S(t), P (t), E(t), I(t), H(t), R(t) in the system (2) remain positive for all t > 0. This completes the proof. Theorem 3.2. (Boundedness) There exists a positively invariant region Ω in which the solution (S(t), P (t), E(t), I(t), H(t), R(t)) of the dy- namical system (2) is bounded. Proof. The positivity has already been estab- lished by Theorem (3.1). For this model the total population is N(t) = S(t)+P (t)+E(t)+ I(t) +H(t) +R(t). Then, we obtain: dN dt = Π− µN − (ρI + ξH). This implies that dN dt ≤ Π− µN, Since the solution I(t) and H(t) are positive. Solving the differential inequality we get the relation, N(t) ≤ Π µ + ( N(0)− Π µ ) e−µt, If N(0) ≤ Π µ , then we obtain 0 ≤ N(t) ≤ Π µ , for all t ≥ 0. If N(0) ≥ Π µ , then we have 0 ≤ N(t) ≤ N(0), for all t ≥ 0. Thus, the feasible solution set of the system (2) remain in the region Ω = { (S, P,E, I,H,R) ∈ R6 + : 0 ≤ N(t) ≤ max ( N(0), Π µ )} . If we start with initial data N(0) ∈ Ω, then the solution N(t) ∈ Ω, for every t > 0. This shows the positively invariance of Ω. Thus, the solution of the dynamical system (2) is bounded. 3.2. Steady state The steady states of the system (2) are so- lutions of the following equations: 0 =Π + ηωR + νP − β(σ1I + σ2E + σ3H)S − (θ + µ)S 0 =θS + (1− η)ωR− (ν + µ)P 0 =β(σ1I + σ2E + σ3H)S − (δ + µ)E 0 =τδE − (α + ε+ µ+ ρ)I 0 =αI − (γ + µ+ ξ)H 0 =(1− τ)δE + εI + γH − (ω + µ)R. There are at most two steady states for the system (2): the disease free equilibrium e0 and endemic equilibrium e1. The disease-free equi- librium point of our model is obtained by set- ting the disease state variables E = 0, I = 0 and H = 0. Thus, the disease free equilibrium point is given by e0 = ( S0, P 0, E0, I0, H0, R0 ) = ( Π(ν + µ) µ(θ + ν + µ) , Πθ µ(θ + ν + µ) , 0, 0, 0, 0 ) . The existence of the endemic equilibrium point depends on the basic reproduction num- ber R0 and will be presented later. 3.3. Basic reproduction number The basic reproduction number, which is denoted by R0, and defined as the average number of secondary infections produced by a single infected individual in a completely sus- East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 63 ceptible population. Using the next genera- tion matrix method Diekmann et al. (2010), the basic reproduction number R0 can be cal- culated from the relation R0 = ρ(FV −1). Let F be the vector for the newly infected and V be the vector for the transfer of individu- als into and out of the infected compartments. Let x = (E, I,H), then we obtain: F(x) =  β(σ1I + σ2E + σ3H)S 0 0  and V(x) =  (δ + µ)E (α + ε+ µ+ ρ)I − τδE (γ + µ+ ξ)H − αI  . The Jacobian matrix to F and V are F = [ ∂Fi(e0) ∂xj ] =  βσ2S 0 βσ1S 0 βσ3S 0 0 0 0 0 0 0  , V = [ ∂Vi(e0) ∂xj ] =  δ + µ 0 0 −τδ k1 0 0 −α k2  , where k1 = α + ε+ µ+ ρ and k2 = γ + µ+ ξ. The next-generation matrix FV −1 is given by FV −1 =  R1 +R2 +R3 βS0(σ1k2+σ3α) k1k2 βσ3S0 k2 0 0 0 0 0 0  , where R1 = βσ2S0 δ + µ , R2 = βσ1S0τδ k1(δ + µ) and R3 = βσ3S0τδα k1k2(δ + µ) . (3) We find the eigenvalues of FV −1 by solving the characteristic equation |FV −1 − λI| = 0 as λ1 = R1 +R2 +R3, λ2 = 0 and λ3 = 0. The basic reproduction number R0 is the spec- tral radius (the largest eigenvalues in modulus) of FV −1 which is given by R0 = ρ(FV −1) = R1 +R2 +R3. The parts R1, R2 and R3 represent the contri- butions from the human-to-human transmis- sion routes (exposed-to-susceptible, infected- to-susceptible and hospitalized-to-susceptible individuals, respectively). We can rewrite the basic reproduction number as follows: R0 = Πβ(ν + µ) µ(θ + ν + µ)(δ + µ) ( σ2 + σ1τδ α + ε+ µ+ ρ + σ3τδα (α + ε+ µ+ ρ)(γ + µ+ ξ) ) . (4) 3.4. Local stability of disease free equilibrium Theorem 3.3. The disease free equilibrium point e0 of the system (2) is locally asymptotically stable if R0 < 1 and unstable if R0 > 1. Proof. The Jacobian matrix of the system (2) at the disease-free equilibrium e0 is given by J(e0) =  −(θ + µ) ν −βσ2S 0 −βσ1S 0 −βσ3S 0 ηω θ −(ν + µ) 0 0 0 (1− η)ω 0 0 βσ2S 0 − (δ + µ) βσ1S 0 βσ3S 0 0 0 0 τδ −k1 0 0 0 0 0 α −k2 0 0 0 (1− τ)δ ε γ −(ω + µ)  . East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 64 The matrix J(e0) is an upper trian- gular block matrix. Its eigenvalues are λ1, λ2, λ3, λ4, λ5 and λ6. Where λ1, λ2 are eigenvalues of the first block matrix of J(e0) and λ3, λ4, λ5, λ6 are eigenvalues of the fourth block matrix of J(e0). The first block matrix of J(e0) given by J1(e0) = [ −(θ + µ) ν θ −(ν + µ) ] . We find the eigenvalues of J1(e0) by solving the characteristic equation |J1(e0)−λI| = 0 as λ1 = −µ and λ2 = −(θ + ν + µ). The fourth block matrix of J(e0) is given by J4(e0) =  βσ2S 0 − (δ + µ) βσ1S 0 βσ3S 0 0 τδ −k1 0 0 0 α −k2 0 (1− τ)δ ε γ −(ω + µ)  . Thus the eigenvalues λ3, λ4, λ5, λ6 are obtained from the characteristic equation of J4(e0): (−(ω + µ)− λ)[λ3 + ((δ + µ)(1−R1) + k1 + k2)λ 2 + [k1(δ + µ)(1− (R1 +R2)) +k2(δ + µ)(1−R1) + k1k2]λ+ k1k2(δ + µ)(1−R0)] = 0. From this equation, we obtain the values for λ to be λ3 = −(ω + µ) and the eigenvalues λ4, λ5, λ6 are the roots of the cubic polynomial: p(λ) = a0λ 3 + a1λ 2 + a2λ+ a3 = 0, where a0 = 1, a1 = (δ + µ)(1−R1) + k1 + k2, a2 = k1(δ + µ)(1− (R1 +R2)) + k2(δ + µ)(1−R1) + k1k2, a3 = k1k2(δ + µ)(1−R0). Furthermore, a1a2 − a3 = k1(δ + µ)2(1−R1)(1− (R1 +R2)) + k2(δ + µ)2(1−R1) 2 + k1k2(δ + µ)(1−R1) +k21(δ + µ)(1−R0) + k1(k1 + k2)[k1(δ + µ)R3 + k2(δ + µ)(1−R1) + k1k2]. If R0 < 1, then R1, R2, R3 and R1 +R2 are strictly less than one. Since R0 = R1+R2+R3. The coefficients a1, a2 and a3 are positive and a1a2 > a3 if R0 < 1. Thus, all the eigen- values of J(e0) are negative. It follows by Routh-Hurwitz criteria that the disease free equilibrium e0 is locally asymptotically stable for R0 < 1. If R0 > 1, then a3 is negative and the Routh-Hurwitz criterion tells that the disease free equilibrium e0 is unstable. East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 65 3.5. Global stability of disease free equilibrium Theorem 3.4. For R0 < 1, the disease free equilibrium e0 of the system (2) is globally asymptot- ically stable if S0 ≥ S. Proof. Let us rewrite our model system (2) as dZ1 dt = F (Z1, Z2), dZ2 dt = G(Z1, Z2), G(Z1, 0) = 0. Where Z1 = (S, P,R) ∈ R3 + represents the class of uninfected individuals and Z2 = (E, I,H) ∈ R3 + represents the class of infected individuals. The disease free equilibrium point of the model is denoted by U0 = (Z∗ 1 , 0), where Z∗ 1 = ( Π(ν+µ) µ(θ+ν+µ) , Πθ µ(θ+ν+µ) , 0 ) . Since the disease free equilibrium point is locally asymptotically stable (see theorem (3.3) ), to prove global sta- bility, we will apply the Castillo-Chavez theo- rem Castillo-Chavez et al. (2002). From sys- tem (2), we have dZ1 dt = F (Z1, Z2) =  Π+ ηωR + νP − β(σ1I + σ2E + σ3H)S − (θ + µ)S θS + (1− η)ωR− (ν + µ)P (1− τ)δE + εI + γH − (ω + µ)R  , dZ2 dt = G(Z1, Z2) =  β(σ1I + σ2E + σ3H)S − (δ + µ)E τδE − (α + ε+ µ+ ρ)I αI − (γ + µ+ ξ)H  . I. To show Z∗ 1 is globally asymptotically stable for the system dZ1 dt = F (Z1, 0), let us consider the reduced system dZ1 dt = F (Z1, 0) =  Π+ ηωR + νP − (θ + µ)S θS + (1− η)ωR− (ν + µ)P −(ω + µ)R  . (5) We can rewrite the system (5) as: dS dt = −(θ + µ)S + νP + ηωR +Π, (61) dP dt = θS − (ν + µ)P + (1− η)ωR, (62) dR dt = −(ω + µ)R, (63) East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 66 and admits as solutions S(t) = Π(ν + µ) µ(θ + ν + µ) + 1 θ + ν [( ν(S(0) + P (0) +R(0))− Πν µ ) e−µt + ( θS(0)− νP (0) + ω(ν − η(θ + ν)) θ + ν − ω R(0)− Πθ θ + ν + µ ) e−(θ+ν+µ)t ] − (ν − ηω) θ + ν − ω R(0)e−(ω+µ)t, P (t) = Πθ µ(θ + ν + µ) + 1 θ + ν [( θ(S(0) + P (0) +R(0))− Πθ µ ) e−µt − ( θS(0)− νP (0) + ω(ν − η(θ + ν)) θ + ν − ω R(0)− Πθ θ + ν + µ ) e−(θ+ν+µ)t ] − (θ − (1− η)ω) θ + ν − ω R(0)e−(ω+µ)t, R(t) = R(0)e−(ω+µ)t. Taking the limit as t goes to ∞, we obtain (S(t), P (t), R(t)) → ( Π(ν + µ) µ(θ + ν + µ) , Πθ µ(θ + ν + µ) , 0 ) = Z∗ 1 . Therefore, Z∗ 1 is globally asymptotically stable for the system dZ1 dt = F (Z1, 0). II. We will show that G(Z1, Z2) = AZ2 − Ĝ(Z1, Z2), Ĝ(Z1, Z2) ≥ 0 for (Z1, Z2) ∈ Ω where A = ∂G ∂Z2 (Z∗ 1 , 0) is a Metzler matrix (the off diagonal elements of A are non-negative) and Ω is the region where the model makes biological sense. Consider a matrix A = ∂G ∂Z2 (Z∗ 1 , 0) =  βσ2S 0 − (δ + µ) βσ1S 0 βσ3S 0 τδ −(α + ε+ µ+ ρ) 0 0 α −(γ + µ+ ξ)  . Hence, A is a Metzer matrix (off diagonal elements are non-negative). Here, Ĝ(Z1, Z2) = AZ2 −G(Z1, Z2). After some simplification, we obtain Ĝ(Z1, Z2) =  βσ2E(S 0 − S) + βσ1I(S 0 − S) + βσ3H(S0 − S) 0 0  , Ĝ(Z1, Z2) = (S0 − S)  β(σ2E + σ1I + σ3H) 0 0  ≥ 0. Therefore by Castillo-Chavez theorem Castillo-Chavez et al. (2002), the disease free equilibrium point e0 of the system (2) is globally asymptotically stable for R0 < 1. East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 67 3.6. Endemic equilibrium point Endemic equilibrium point is a steady state solution where the disease persists in the pop- ulation. In the presence of disease in the population, there exist an equilibrium point called endemic equilibrium point denoted by e1 = (S∗, P ∗, E∗, I∗, H∗, R∗). It can be ob- tained by setting each equation of the system (2) equal to zero. Then we obtained S∗ = S0 R0 , P ∗ = θS0K + (1− η)ωµδ(θ + ν + µ)[(1− τ)k1k2 + ετk2 + γατ ]S0(R0 − 1) (ν + µ)KR0 , E∗ = k1k2µ(ω + µ)(θ + ν + µ)S0(R0 − 1) KR0 , I∗ = k2δµτ(ω + µ)(θ + ν + µ)S0(R0 − 1) KR0 , H∗ = αδµτ(ω + µ)(θ + ν + µ)S0(R0 − 1) KR0 , R∗ = δµ(θ + ν + µ)[(1− τ)k1k2 + ετk2 + γατ ]S0(R0 − 1) KR0 , where K = k1k2[µ(ν + µ)(ω + δ + µ) + δµ(1− η)ω] + ωδτ(ηµ+ ν)[α(µ+ ξ) + k2(µ+ ρ)], provided that R0 > 1. From this we see that for the endemic equilibrium to exist R0 > 1. Moreover, the force of infection can be updated as λ∗ = β(σ1I ∗ + σ2E ∗ + σ3H ∗). (7) When we substitute the expression for E∗, I∗ and H∗ into the force of infection λ∗, we obtain λ∗ = µk1k2(δ + µ)(ω + µ)(θ + ν + µ)(R0 − 1) K , (8) provided that R0 > 1. From this, we see that, there is no endemic equilibrium of the system (2) if R0 < 1. Therefore, this condi- tion shows that it is not possible for backward bifurcation in this model. Hence we have es- tablished the following result. 3.7. Local stability of endemic equilib- rium Theorem 3.6. The endemic equilibrium e1 of the system (2) is locally asymptotically stable Theorem 3.5. A unique endemic equilibrium point e1 = (S∗, P ∗, E∗, I∗, H∗, R∗) exists and is positive if R0 > 1. When we plot the force of infection λ∗ over R0 by using the expression for λ∗ we got a for- ward bifurcation in Figure (2). if R0 > 1. East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 68 Proof. To determine the local stability of en- demic equilibrium, we used the center manifold theory Castillo-Chavez and Song (2004), by taking β as a bifurcation parameter. We make the following change of variables on the system (2). Let S = x1, P = x2, E = x3, I = x4, H = x5 and R = x6. Moreover, by using vector no- tation x = (x1, x2, x3, x4, x5, x6) T , the system (2) can be written in the form dx dt = F (x), with F = (f1, f2, f3, f4, f5, f6) T . We choose β = β∗ as a bifurcation param- eter. Solving for β∗ from R0 = 1, we obtain β∗ = µ(θ + ν + µ)(δ + µ)(α + ε+ µ+ ρ)(γ + µ+ ξ) Π(ν + µ)[σ2(α + ε+ µ+ ρ)(γ + µ+ ξ) + σ1τδ(γ + µ+ ξ) + σ3τδα] . The Jacobian matrix of the system (2) evaluated at the disease free equilibrium e0 with β = β∗ is given by J∗ =  −(θ + µ) ν −β∗σ2S 0 −β∗σ1S 0 −β∗σ3S 0 ηω θ −(ν + µ) 0 0 0 (1− η)ω 0 0 β∗σ2S 0 − (δ + µ) β∗σ1S 0 β∗σ3S 0 0 0 0 τδ −k1 0 0 0 0 0 α −k2 0 0 0 (1− τ)δ ε γ −(ω + µ)  , where k1 = α + ε+ µ+ ρ, k2 = γ + µ+ ξ. The Jacobian matrix J∗ of the linearized system has a simple zero eigenvalue with all other eigenvalues having negative real part, hence the center manifold theory will be used to analyse the dynamics of the sys- tem near β = β∗. Thus, e0 is a non- hyperbolic equilibrium, when β = β∗. Now, the components of the right eigenvector w = (w1, w2, w3, w4, w5, w6) T of J∗ associated with the zero eigenvalue are given by w1 = −k1k2[µ(ν + µ)(ω + δ + µ) + δµ(1− η)ω] + ωδτ(ηµ+ ν)[α(µ+ ξ) + k2(µ+ ρ)] τδµk2(θ + ν + µ)(ω + µ) w4, w2 = −ωδµ[τ (k2(µ+ ρ) + α(µ+ ξ)) (2− η)− (1− η)] + θµ(δ + ω + µ) τδµk2(θ + ν + µ)(ω + µ) w4, w3 = α + ε+ µ+ ρ τδ w4, w4 = w4 > 0, w5 = α γ + µ+ ξ w4, w6 = (1− τ)(α + ε+ µ+ ρ)(γ + µ+ ξ) + ετ(γ + µ+ ξ) + γατ τ(ω + µ)(γ + µ+ ξ) w4. Similarly, the components of the left eigen- vector v = (v1, v2, v3, v4, v5, v6) T of J∗ associ- ated with the zero eigenvalue are given by v1 = v2 = v6 = 0, v3 = v3 > 0, v4 = β∗S0(σ1(γ + µ+ ξ) + ασ3) (α + ε+ µ+ ρ)(γ + µ+ ξ) v3, v5 = β∗S0σ3 γ + µ+ ξ v3. East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 69 Since the first, second and six component of v are zero, we don’t need the partial derivatives of f1, f2 and f6. From the partial derivatives of f3, f4 and f5 at the disease free equilibrium point, the only ones that are nonzero are the following: ∂2f3 ∂x3∂x1 = ∂2f3 ∂x1∂x3 = β∗σ2, ∂2f3 ∂x4∂x1 = ∂2f3 ∂x1∂x4 = β∗σ1, ∂2f3 ∂x5∂x1 = ∂2f3 ∂x1∂x5 = β∗σ3, ∂2f3 ∂x3∂β = σ2S 0, ∂2f3 ∂x4∂β = σ1S 0, ∂2f3 ∂x5∂β = σ3S 0. The direction of the bifurcation at R0 = 1 is determined by the signs of the bifurcation coefficients a and b. Hence, a = v3 6∑ i,j=1 wiwj ∂2f3 ∂xi∂xj (S0, P 0, 0, 0, 0, 0) = 2v3w1 (w3β ∗σ2 + w4β ∗σ1 + w5β ∗σ3) = 2β∗w1 ( σ2k1 τδ + σ1 + ασ3 k2 ) v3w4 =− 2(δ + µ) (τδ)2k2Π(ν + µ)(ω + µ) [ k21k2 (µ(ν + µ)(ω + δ + µ) + δµ(1− η)ω) + ωδτk1(ηµ+ ν) (α(µ+ ξ) + k2(µ+ ρ))] v3w 2 4 < 0. and b = v3 6∑ i=1 wi ∂2f3 ∂xi∂β (S0, P 0, 0, 0, 0, 0) = v3 ( k1w4σ2S 0 τδ + w4σ1S 0 + αw4σ3S 0 k2 ) = (δ + µ)(α + ε+ µ+ ρ) β∗τδ v3w4 > 0. Since a < 0 and b > 0 at β = β∗. Based on the Theorem 4.1 stated in Castillo-Chavez and Song (2004), the system (2) undergoes a for- ward bifurcation at R0 = 1 and the unique en- demic equilibrium e1 is locally asymptotically stable for R0 > 1. 3.8. Bifurcation analysis We investigate the nature of the bifurcation by using the center manifold theory Castillo- Chavez and Song (2004). In short, the the- ory is summarized by Theorem 4.1 in Castillo- Chavez and Song (2004). In such a theorem, there are two important quantities: the coef- ficients, say a and b, of the normal form rep- resenting the dynamics of the system on the central manifold. These coefficients decide the bifurcation. In particular, if a < 0 and b > 0, then the bifurcation is forward. In the proof of Theorem (3.6), we have already justified the system (2) undergoes a forward bifurcation at R0 = 1. Thus the basic reproduction num- ber R0 plays an important role in the disease spread. If R0 < 1, then its easy to control the disease but if R0 > 1, then the society will ex- perience endemic disease spreading. The for- ward bifurcation diagram can be seen in Figure (2). East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 70 Figure 2: Forward bifurcation diagram for the COVID-19 model (2). From Figure (2), it is clear that when R0 < 1, the system (2) has no endemic equilib- rium and the disease-free equilibrium is stable. When R0 > 1, a stable endemic equilibrium appears and the disease free equilibrium be- comes unstable, i.e. exchange of stability of the equilibrium’s (forward bifurcation) occurs at the bifurcation point R∗ 0 = 1. 3.9. Sensitivity analysis Sensitivity analysis is a useful tool in model building as well as in model evaluation by showing how the model behavior responds to changes in parameter values Martcheva (2015). The threshold parameter R0 which de- termines stability is a function of the param- eters Π, β, σ1, σ2, σ3, ν, µ, θ, δ, γ, ξ, α, ε, ρ, τ. We recall that the basic reproduction number R0 is given by R0 = Πβ(ν + µ) [(γ + µ+ ξ) [σ2(α + ε+ µ+ ρ) + σ1τδ] + σ3τδα] µ(θ + ν + µ)(δ + µ)(α + ε+ µ+ ρ)(γ + µ+ ξ) . Thus, in order to identify the most sensi- tive parameters for model (2), we compute the relative sensitivity of R0 with respect to the above parameters. Then using the parameter values from Table (3), we display the sensitiv- ity indices of R0 with respect to the parameters in Figure (3). East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 71 △R0 β = ∂R0 ∂β × β R0 = 1, △R0 Π = ∂R0 ∂Π × Π R0 = 1, △R0 σ1 = ∂R0 ∂σ1 × σ1 R0 = σ1τδ(γ + µ+ ξ) (γ + µ+ ξ) [σ2(α + ε+ µ+ ρ) + σ1τδ] + σ3τδα , △R0 σ2 = ∂R0 ∂σ2 × σ2 R0 = σ2(α + ε+ µ+ ρ)(γ + µ+ ξ) (γ + µ+ ξ) [σ2(α + ε+ µ+ ρ) + σ1τδ] + σ3τδα , △R0 σ3 = ∂R0 ∂σ3 × σ3 R0 = σ3τδα (γ + µ+ ξ) [σ2(α + ε+ µ+ ρ) + σ1τδ] + σ3τδα , △R0 τ = ∂R0 ∂τ × τ R0 = τδ [σ1(γ + µ+ ξ) + σ3α] (γ + µ+ ξ) [σ2(α + ε+ µ+ ρ) + σ1τδ] + σ3τδα , △R0 ν = ∂R0 ∂ν × ν R0 = θν (θ + ν + µ)(ν + µ) , △R0 α = ∂R0 ∂α × α R0 = ατδ [σ3(ε+ µ+ ρ)− σ1(γ + µ+ ξ)] (α + ε+ µ+ ρ) [(γ + µ+ ξ) [σ2(α + ε+ µ+ ρ) + σ1τδ] + σ3τδα] , △R0 δ = ∂R0 ∂δ × δ R0 = δ [(γ + µ+ ξ) [σ1τµ− σ2(α + ε+ µ+ ρ)] + σ3ταµ] (δ + µ) [(γ + µ+ ξ) [σ2(α + ε+ µ+ ρ) + σ1τδ] + σ3τδα] , △R0 θ = ∂R0 ∂θ × θ R0 = − θ θ + ν + µ , △R0 ε = ∂R0 ∂ε × ε R0 = − ετδ [σ1(γ + µ+ ξ) + σ3α] (α + ε+ µ+ ρ) [(γ + µ+ ξ) [σ2(α + ε+ µ+ ρ) + σ1τδ] + σ3τδα] , △R0 ρ = ∂R0 ∂ρ × ρ R0 = − ρτδ [σ1(γ + µ+ ξ) + σ3α] (α + ε+ µ+ ρ) [(γ + µ+ ξ) [σ2(α + ε+ µ+ ρ) + σ1τδ] + σ3τδα] , △R0 γ = ∂R0 ∂γ × γ R0 = − σ3τδαγ (γ + µ+ ξ) [(γ + µ+ ξ) [σ2(α + ε+ µ+ ρ) + σ1τδ] + σ3τδα] , △R0 ξ = ∂R0 ∂ξ × ξ R0 = − σ3τδαξ (γ + µ+ ξ) [(γ + µ+ ξ) [σ2(α + ε+ µ+ ρ) + σ1τδ] + σ3τδα] , △R0 µ = ∂R0 ∂µ × µ R0 = 1 (ν + µ)D [ µD − [ σ2A(α + ε+ µ+ ρ)2(γ + µ+ ξ)2(ν + µ) + σ1τδ[(α + ε+ ρ)A+B](γ + µ+ ξ)2(ν + µ) + σ3τδα [(γ + ξ)(A+B) + (α + ε+ ρ)B + C] (ν + µ)]] , where A = δ(θ + ν) + 2µ(δ + θ + ν) + 3µ2, B = 2µδ(θ + ν) + 3µ2(δ + θ + ν) + 4µ3, C = 3µ2δ(θ + ν) + 4µ(δ + θ + ν) + 5µ4, D = (α + ε+ µ+ ρ)(γ + µ+ ξ)(δ + µ)(θ + ν + µ) [(γ + µ+ ξ) [σ2(α + ε+ µ+ ρ) + σ1τδ] + σ3τδα] . East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 72 Note that the sensitivity index may depend on several parameters of the system, but also can be constant, independent of any parame- ter. For example, △R0 β = +1 means that in- creasing (decreasing) β by a given percentage increases (decreases) always R0 by that same percentage. Figure 3: The sensitivity indices of R0 with respect to the parameters. Figure (3) shows that the recruitment rate Π, the contact rate β, the modification param- eter σ2, the exposed progression rate δ, the natural death rate µ , the protection rate θ and the waning rate of protected individuals to sus- ceptible class ν are the most sensitive param- eters for R0. The parameters τ, σ1, σ2, σ3, ν, β and Π have positive correlation with R0. This indicates that the spread of COVID-19 de- creases with decrease of these parameters. The parameters ξ, θ, ρ, µ, γ, ε, δ and α have nega- tive correlation with R0. This implies that the spread of the virus decreases with an increase of these parameters. 4. NUMERICAL SIMULATIONS AND DISCUSSION In this section, we perform numerical sim- ulation to support our analytical results. The numerical simulations are carried out with help of the ode45 Matlab tool. Using the param- eter values given in Table (3) and the ini- tial conditions below in the model equations (2) simulation study is conducted. The pa- rameter values have been taken based on the literature and the real characteristics of the virus. Furthermore, initial conditions are de- termined as follow. The total population of Ethiopia for the year 2021 is estimated about N(0) = 114, 963, 588 people Kifle and Obsu (2022). On December 31, 2020 the total ac- tive cases (Infected individuals are 10,245 i.e I(0) = 10, 245), and the Hospitalized are H(0) = 5062 and the total recovered by the date are R(0) = 81, 144 and we are assumed P (0) = 108076, E(0) = 15000 and hence is S(0) = N(0)− (P (0) + E(0) + I(0) +H(0) + R(0)) = 114, 744, 061. East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 73 Table 3: The parameter values of the modified model (per day). Parameter Value Source Π 1300 Assumed θ 0.7 Assumed ν 0.15 Assumed β 0.00000058 Assumed σ1 0.0001 Alemneh and Telahun (2020) σ2 0.02 Alemneh and Telahun (2020) σ3 0.00003 Assumed δ 1/14 Kifle and Obsu (2022) α 0.04 Assumed µ 1/(64*12*30)Kifle and Obsu (2022) ρ 0.0004 Alemneh and Telahun (2020) ξ 0.015 Assumed τ 0.7 Alemneh and Telahun (2020) ε 0.0476 Assumed γ 0.033 Assumed η 0.4 Assumed ω 0.011 Assumed East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 74 Figure (4) shows the predicted total con- firmed cases and the real data total confirmed cases for Ethiopia from 31 December 2020 to 30 March 2022 (It the data has been taken COVID-19 pandemic (2022) ). There is some differences between the prediction and the real data that should come due to the fact that there were no enough covid test kits during the early time not only in Ethiopia but also all over the world. After wards we continue to validate the local stability of DFE and EEP of model. Figure 4: The prediction using model (2) and the real data of confirmed cases for Ethiopia from 31 December 2020 to 30 March 2022 In Figure (5) with R0 = 0.22950, we ob- serve that for the basic reproduction number R0 < 1, all solutions curve goes to the disease free equilibrium point. As a result, the disease goes to extinct or the disease dies out. In Figure (6) with R0 = 1.10791, we ob- serve that for the basic reproduction number R0 > 1, all solutions curve goes away from the disease free equilibrium point. These in- dicate that the disease-free equilibrium point is unstable for the values of R0 > 1, and the solutions will go to the endemic equilibrium point. Consequently, the disease invade in a population. In Figure (7), we observe as the protection rate θ increases, all infected classes will signif- icantly decrease over time. This confirmed the result from the fact that strict use of safety (protection) measures within the population plays a critical role in confine the spread of the disease as in Bachar et al. (2021). It is pre- dicted that the population will be disease-free. This is due to the case that when protection rate increases the basic reproduction number decreases. East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 75 Figure 5: The time series plot of model (2) when R0 = 0.22950 In Figure (8), we observe that as the con- tact rate β decrease, R0 and also all infected classes are decreasing. Further this habitual the result obtained from the fact that a de- crease in contact among the population plays a indispensable role in curtailing the spread of the disease as in Ahmed et al. (2021). 5. EXTENSION OF THE MODIFIED MODEL INTO AN OPTIMAL CON- TROL In this section, we will use optimal con- trol theory to find protection (for suscepti- ble) and hospitalization (for infected) strate- gies that would mitigate the spread of COVID- 19 in the population. 5.1. Optimal protection and hospitaliza- tion using modified model In this subsection, we study the optimal protection of susceptible population and hos- pitalization of infected individuals in order to minimize the outbreak of COVID-19 in the population. Let us define our control set U to be U = {(θ(t), α(t)) : 0 ≤ θ(t), α(t) ≤ ϵi, 0 ≤ t ≤ T, 0 < ϵi ≤ 1, i = 1, 2}, (9) where θ(t) and α(t) are Lebesgue measurable quantities bounded above by ϵ1 and ϵ2 respec- tively. We will minimize the objective func- tional J [θ(t), α(t)] = ∫ T 0 ( B1E(t) +B2I(t) + 1 2 ( A1θ 2(t) +A2α 2(t) )) dt, (10) where constants B1, B2, A1 and A2 are posi- tive. Here, we want to find the optimal val- ues θ(t) and α(t) that minimizes the objective functional (10) subject to the state system (2). The goal is to find the optimal control (θ∗(t), α∗(t)) such that J [θ∗(t), α∗(t)] = min (θ,α)∈U J [θ(t), α(t)]. (11) East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 76 Figure 6: The time series plot of model (2) when β = 0.0000028 with R0 = 1.11073 Figure 7: Impact of protection rate θ on R0 and infected classes. 5.2. Existence of an optimal control The existence of the optimal control can be showed by using an approach of Fleming and East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 77 Figure 8: Impact of contact rate β on R0 and infected classes. Rishel (2012). Theorem 5.1. Given the objective functional J(θ(t), α(t)) (10) with admissible control set U , subject to the state system (2), then there exist an optimal control double u∗ = (θ∗, α∗) in U such that J [θ∗, α∗] = min (θ,α)∈U J [θ(t), α(t)]. (12) Proof. To prove the existence of optimal con- trol, we need to verify the following conditions. (a) The set of solutions to the state system (2) and control parameters in (9) are non-empty. (b) The set U is convex and closed. (c) The right hand side of system (2) is bounded above by sum of bounded con- trol and state and can be written as a lin- ear function of the control variables with coefficients dependent on time and state variables. (d) The integrand function L(E, I, θ, α, t) is convex on U and L(E, I, θ, α, t) ≥ h(u), where h(u) is continuous and ||u||−1h(u) → ∞ when ||u|| → ∞. Here u = (θ(t), α(t)). In Theorem (3.2), we have already justified the boundedness of the solution of the state sys- tem (2). Since our solution for the model is bounded by N(t) ≤ max ( N(0), Π µ ) for all t ≥ 0. This implies that the solutions of the state system are continuous and bounded for each admissible control functions in U . Moreover, the right hand side of the model equations (2) satisfes the Lipschitz condition with respect to state variables. Hence, the state system (2) has a unique solution corresponding to each admissible control function (θ(t), α(t)) ∈ U . Thus, condition (a) is achieved. To verify condition (b), given that the con- trol set U = {u ∈ R2 : ||u||∞ ≤ 1}. Let ψ ∈ [0, 1] and v1, v2 ∈ U such that ||v1||∞ ≤ 1 and ||v2||∞ ≤ 1, then ||ψv1+(1−ψ)v2||∞ ≤ ψ||v1||∞+(1−ψ)||v2||∞ ≤ 1. Thus, the set U is convex and closed. East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 78 To verify condition (c), let u = (θ(t), α(t)) ∈ U, X = (S, P,E, I,H,R) and the right hand side of the state system (2) is given by f(t,X, u) =  Π+ ηωR + νP − β(σ1I + σ2E + σ3H)S − (θ(t) + µ)S θ(t)S + (1− η)ωR− (ν + µ)P β(σ1I + σ2E + σ3H)S − (δ + µ)E τδE − (α(t) + ε+ µ+ ρ)I α(t)I − (γ + µ+ ξ)H (1− τ)δE + εI + γH − (ω + µ)R  . (13) Then from (13) we get, f(t,X, u) = g(t,X) + h(t,X)uT , where g(t,X) =  Π+ ηωR + νP − β(σ1I + σ2E + σ3H)S − µS (1− η)ωR− (ν + µ)P β(σ1I + σ2E + σ3H)S − (δ + µ)E τδE − (ε+ µ+ ρ)I −(γ + µ+ ξ)H (1− τ)δE + εI + γH − (ω + µ)R  and h(t,X) =  −S 0 S 0 0 0 0 −I 0 I 0 0  . Since, by using the properties of a norm of a matrix we have, ||f(t,X, u)|| = ||g(t,X) + h(t,X)uT || ≤ ||g(t,X)||+ ||h(t,X)|| ||u||. Thus, condition (c) is proved. To verify condition (d), the integrand of the objective functional (10) L(E, I, θ, α, t) = B1E(t) +B2I(t) + 1 2 ( A1θ 2(t) + A2α 2(t) ) (14) is the sum of convex function and hence con- vex with respect to control parameters θ(t) and α(t). Moreover, L(E, I, θ, α, t) = B1E(t) +B2I(t) + 1 2 ( A1θ 2(t) + A2α 2(t) ) ≥ 1 2 ( A1θ 2(t) + A2α 2(t) ) . We define a continuous function h(u) = ϕ||u||2, where ϕ = min ( A1 2 , A2 2 ) > 0 and u = (θ(t), α(t)). Then we have L(E, I, θ, α, t) ≥ 1 2 ( A1θ 2(t) +A2α 2(t) ) ≥ ϕ||u||2, (15) since ϕ = min ( A1 2 , A2 2 ) > 0. This im- plies that L(E, I, θ, α, t) ≥ h(u). Consider, ||u||−1h(u) = ||u||−1ϕ||u||2 = ϕ||u||. This gives that ||u||−1h(u) = ϕ||u|| → ∞ when ||u|| → ∞. Thus, condition (d) is proved. Hence, all conditions (a)-(d) shows that there exists an optimal control u∗ = (θ∗, α∗) that minimizes the cost functional J(θ(t), α(t)) over U . There- fore, the existence of optimal control is estab- lished. 5.3. The Hamiltonian and optimality system We used Pontryangin’s Maximum Princi- ple Lenhart and Workman (2007) to drive the necessary conditions that an optimal control must satisfy. This principle converts the objec- tive functional (10) subject to the state system (2) into a problem of minimizing point-wise a Hamiltonian (H), with respect to θ(t) and α(t) as: East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 79 H = B1E +B2I + 1 2 A1θ 2 + 1 2 A2α 2 + λ1[Π + ηωR + νP − β(σ2E + σ1I + σ3H)S − (θ + µ)S] + λ2[θS + (1− η)ωR− (ν + µ)P ] + λ3[β(σ2E + σ1I + σ3H)S − (δ + µ)E] + λ4[τδE − (α + ε+ µ+ ρ)I] + λ5[αI − (γ + µ+ ξ)H] + λ6[(1− τ)δE + εI + γH − (ω + µ)R], (16) where λi, i = 1, 2, 3, 4, 5, 6, represent the adjoint variables associated with the state variables S, P,E, I,H and R to be determined suitably by applying Pontryagin’s Maximal Principle Lenhart and Workman (2007). Theorem 5.2. For an optimal control set θ, α that minimizes J over U , there are adjoint variables, λ1, ..., λ6 such that: dλ1 dt = [β(σ2E + σ1I + σ3H) + θ + µ]λ1 − θλ2 − β(σ2E + σ1I + σ3H)λ3, dλ2 dt = −νλ1 + (ν + µ)λ2, dλ3 dt = −B1 + βσ2Sλ1 − [βσ2S − (δ + µ)]λ3 − τδλ4 − (1− τ)δλ6, dλ4 dt = −B2 + βσ1Sλ1 − βσ1Sλ3 + (α + ε+ µ+ ρ)λ4 − αλ5 − ελ6, dλ5 dt = βσ3Sλ1 − βσ3Sλ3 + (γ + µ+ ξ)λ5 − γλ6, dλ6 dt = −ηωλ1 − (1− η)ωλ2 + (ω + µ)λ6. (17) with transiversality conditions λi(T ) = 0, i = 1, ..., 6. (18) Moreover, we obtain the control set (θ∗, α∗) characterized by θ∗(t) = max { 0,min ( ϵ1, S(λ1 − λ2) A1 )} , α∗(t) = max { 0,min ( ϵ2, I(λ4 − λ5) A2 )} . (19) Proof. The form of the adjoint equations and transversality conditions are standard results from Pontryagin’s Maximum Principle Lenhart and Workman (2007). We differentiate Hamiltonian (H) (16) with respect to the state variables S, P,E, I,H and R, respectively, and then the adjoint East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 80 system can be written as dλ1 dt = −∂H ∂S = [β(σ2E + σ1I + σ3H) + θ + µ]λ1 − θλ2 − β(σ2E + σ1I + σ3H)λ3, dλ2 dt = −∂H ∂P = −νλ1 + (ν + µ)λ2, dλ3 dt = −∂H ∂E = −B1 + βσ2Sλ1 − [βσ2S − (δ + µ)]λ3 − τδλ4 − (1− τ)δλ6, dλ4 dt = −∂H ∂I = −B2 + βσ1Sλ1 − βσ1Sλ3 + (α + ε+ µ+ ρ)λ4 − αλ5 − ελ6, dλ5 dt = −∂H ∂H = βσ3Sλ1 − βσ3Sλ3 + (γ + µ+ ξ)λ5 − γλ6, dλ6 dt = −∂H ∂R = −ηωλ1 − (1− η)ωλ2 + (ω + µ)λ6. with transversality conditions λi(T ) = 0, i = 1, ..., 6. Similarly by following the approach of Pontryagin et al Pontryagin (2018), the characterization of optimal controls θ∗(t), α∗(t), that is, the optimality equations are obtained based on the conditions: ∂H ∂θ = 0 and ∂H ∂α = 0, which gives, θ = S(λ1 − λ2) A1 , α = I(λ4 − λ5) A2 . Since θ and α are bounded in U by ϵ1 and ϵ2 respectively. Therefore, the optimal controls θ∗(t) and α∗(t) are given by θ∗(t) = max { 0,min ( ϵ1, S(λ1 − λ2) A1 )} , α∗(t) = max { 0,min ( ϵ2, I(λ4 − λ5) A2 )} . This completes the proof. The optimality system is formed from the state system (2) and the adjoint variable sys- tem (17) by incorporating the characterized control set and initial and transversal condi- tion. Then we have the following optimality system: East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 81 dS dt = Π+ ηωR + νP − β(σ1I + σ2E + σ3H)S − (θ∗ + µ)S, dP dt = θ∗S + (1− η)ωR− (ν + µ)P, dE dt = β(σ1I + σ2E + σ3H)S − (δ + µ)E, dI dt = τδE − (α∗ + ε+ µ+ ρ)I, dH dt = α∗I − (γ + µ+ ξ)H, dR dt = (1− τ)δE + εI + γH − (ω + µ)R, dλ1 dt = [β(σ2E + σ1I + σ3H) + θ∗ + µ]λ1 − θ∗λ2 − β(σ2E + σ1I + σ3H)λ3, dλ2 dt = −νλ1 + (ν + µ)λ2, dλ3 dt = −B1 + βσ2Sλ1 − [βσ2S − (δ + µ)]λ3 − τδλ4 − (1− τ)δλ6, dλ4 dt = −B2 + βσ1Sλ1 − βσ1Sλ3 + (α∗ + ε+ µ+ ρ)λ4 − α∗λ5 − ελ6, dλ5 dt = βσ3Sλ1 − βσ3Sλ3 + (γ + µ+ ξ)λ5 − γλ6, dλ6 dt = −ηωλ1 − (1− η)ωλ2 + (ω + µ)λ6, λi(T ) = 0, i = 1, ..., 6, (S(0), P (0), E(0), I(0), H(0), R(0)) = (S0, P0, E0, I0, H0, R0) . 6. Numerical simulations of optimal control problem In this section, we perform some numeri- cal solutions on the modified model (2) and the resulting optimality system consisting of the state equations (2) and the adjoint system (17) with the characterizations (19). We make use of the parameter values given in Table (3) for the simulation. An iterative scheme is used to find the op- timal solution of the optimality system. Since the state system (2) has initial conditions and the adjoint systems (17) have final conditions, we solve the state system using a forward fourth-order Runge–Kutta method and solve the adjoint system using a backward fourth- order Runge–Kutta method ?. The solution iterative scheme involves making a guess of the controls and using that guess to solve the state system. The initial guess of the control to- gether with the solution of the state systems is used to solve the adjoint systems. The con- trols are then updated using a convex combi- nation of the previous controls and the values obtained using the characterizations. The up- dated controls are then used to repeat the so- lution of the state and adjoint systems. This process is repeated until the values in the cur- rent iteration are close enough to the previous iteration values Lenhart and Workman (2007). We used B1 = B2 = 1, A1 = 40, A2 = 80 and final intervention time T = 350 days for simulation of COVID-19 model with op- timal control. Additionally, we used S(0) = 110079, P (0) = 108076, E(0) = 15000, I(0) = 13813, H(0) = 10050, R(0) = 12156 as initial values. East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 82 6.1. Optimal control comparisons and strategies In this subsection, we compare the results of constant and optimal control as did in ?. We first compare the cost of infection for each strategies and then compare the exposed and infected population. The optimal control cost for each parame- ter is less than the constant control at all times as depicted in Figures (91) and (101). It can be observed that for protection rate θ the cost reduction is significant compared to the other parameters which is due to the high sensitivity of parameter over the system. 0 50 100 150 200 250 300 350 time in days 2 4 6 8 10 12 14 16 Co st of Di se as e C on tro l #104 Optimal Control 3=0.01 3=0.24 3=0.47 3=0.7 (1) 0 50 100 150 200 250 300 350 time in days 0 2 4 6 8 10 12 14 Ex po se d P op ula tio n #104 Optimal Control 3=0.01 3=0.24 3=0.47 3=0.7 (2) 0 50 100 150 200 250 300 350 time in days 1 2 3 4 5 6 7 Inf ec ted P op ula tio n #104 Optimal Control 3=0.01 3=0.24 3=0.47 3=0.7 (3) 0 50 100 150 200 250 300 350 time in days 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Pr op or tio n o f P op ula tio n P ro tec ted -= 0.0001 -=0.00015 -= 0.0002 -=0.00025 (d) Figure 9: Optimal control for θ. (a) Cost comparison for optimal and constant controls for θ. (b) Exposed population comparison for optimal and constant controls for θ. (c) Infected population comparison for optimal and constant controls for θ. (d) Optimal strategies for θ with different β. Now we compare the exposed population using constant and optimal controls. In Figure (92), we observe that the optimal time depen- dent protection strategy, θ, yields a significant drop in exposed population when compared to its constant counterparts, similarly using the optimal time dependent hospitalization strat- egy α we see a significant drop in the exposed population as compared to constant rates of hospitalization over time as depicted in Figure (102). We compare the infected population using constant and optimal controls. In Figure (93), we observe that the optimal time dependent protection strategy, θ, yields a significant drop in infected population when compared to its constant counterparts, similarly using the opti- mal time dependent hospitalization strategy α we see a significant drop in the infected popula- tion as compared to constant rates of hospital- ization over time as depicted in Figure (103). As contact rate β increases, so does the time for which the maximum control should be applied. In Figures (9d) and (10d), we ob- East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 83 0 50 100 150 200 250 300 350 time in days 2 3 4 5 6 7 8 9 10 Co st of Di se as e C on tro l #104 Optimal Control ,=0.1 ,=0.16667 ,=0.23333 ,=0.3 (1) 0 50 100 150 200 250 300 350 time in days 1 2 3 4 5 6 7 Ex po se d P op ula tio n #104 Optimal Control ,=0.1 ,=0.16667 ,=0.23333 ,=0.3 (2) 0 50 100 150 200 250 300 350 time in days 0 0.5 1 1.5 2 2.5 3 Inf ec ted P op ula tio n #104 Optimal Control ,=0.1 ,=0.16667 ,=0.23333 ,=0.3 (3) 0 50 100 150 200 250 300 350 time in days 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Pr op or tio n o f P op ula tio n H os pit ali ze d -= 0.0001 -=0.00015 -= 0.0002 -=0.00025 (d) Figure 10: Optimal control for α. (a) Cost comparison for optimal and constant controls for α. (b) Exposed population comparison for optimal and constant controls for α. (c) Infected population comparison for optimal and constant controls for α. (d) Optimal strategies for α with different β. serve that for both protection and hospitaliza- tion the maximum possible protection and hos- pitalization rates should be maintained to the first few days of the disease which can then be reduced over time. Finally, from Figures (112) and (113), one can easily conclude that the combination of the two controls is significantly more effective in reducing the spread of the virus than when each control is singly applied. Hence, the best choice is to apply two controls all together to mitigate the spread of COVID-19 in the population. From Figure (111), one can ob- serve that the combined implementation of the two control measures is the most cost-effective when compared with the single implementa- tion of each control measure. 7. CONCLUSION In this paper, we proposed a deterministic compartmental model to study the transmis- sion dynamics of COVID-19. The model was an extension of the existing SEIR model by in- cluding protected and hospitalized individuals. We established the well-posedness of the modi- fied model by proving the existence, positivity, and boundedness of the solutions. We computed the steady states and the ba- sic reproduction number R0. Based on the reproduction number R0, it is revealed that whenever R0 < 1, the system has only disease free equilibrium e0 which is locally as well as globally asymptotically stable. When R0 > 1, the system has a unique endemic equilibrium e1 which is locally stable and the disease free equilibrium e0 becomes unstable. We have ob- served that the outbreak of the disease dies out if R0 < 1 and the disease is endemic if R0 > 1. Using center manifold theory, bifurcation anal- ysis of the modified model was proven and the model exhibits forward bifurcation at R0 = 1. In addition from sensitivity analysis of R0, we observed that the recruitment rate Π and East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 84 0 50 100 150 200 250 300 350 time in days 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 Co st of Di se as e C on tro l #104 Both Optimal Only , Optimal Only 3 Optimal (1) 0 50 100 150 200 250 300 350 time in days 1 2 3 4 5 6 7 Ex po se d P op ula tio n #104 Both Optimal Only , Optimal Only 3 Optimal (2) 0 50 100 150 200 250 300 350 time in days 0 0.5 1 1.5 2 2.5 3 3.5 Inf ec ted P op ula tio n #104 Both Optimal Only , Optimal Only 3 Optimal (3) Figure 11: Optimal control for θ, α and both. (a) Cost comparison for optimal θ, α (b) Exposedand both. population comparison. (c) Infected population comparison. References Ahmed I., Modu G. U., Yusuf A., Ku- mam P. and Yusuf I. (2021), ‘A math- ematical model of coronavirus disease (covid-19) containing asymptomatic and symptomatic classes’, Results in physics 21, 103776. Alemneh H. T. and Telahun G. T. (2020), ‘Mathematical modeling and optimal con- trol analysis of covid-19 in ethiopia’, medRxiv . Allen L., Brauer F., van den Driessche P. and Wu, J. (2008), ‘Mathematical epidemiol- ogy, volume 1945 of lecture notes in math- ematics’, Springer, Berlin 13, 14. Bachar M., Khamsi M. A. and Bounkhel M. (2021), ‘A mathematical model for the spread of covid-19 and control mech- anisms in saudi arabia’, Advances in Dif- ference Equations 2021(1), 1–18. Brauer F., Castillo-Chavez C. and Feng Z. East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 85 contact rate β are most sensitive parameters to our model. Numerical results support the fact that decrease in the contact rate β causes the decrease in the value of R0 and after a certain level of β, R0 become less than one. If the protection rate increases, then all in- fected classes are decrease. It is predicted that the population will be disease-free. Furthermore, using optimal control the- ory we suggest protection and hospitaliza- tion strategies. Pontryagin’s Maximum Principle is used to establish the existence and characterization of optimal controls. The study demonstrates that the combined application of both controls is much more effective in reducing the spread of the virus than when each control is applied individually. Finally, our analysis indicates that an optimal control is more preferable than maintaining a high constant control. (2019), Mathematical models in epidemi- ology, Vol. 32, Springer. Castillo-Chavez C. Blower, S. Van den Driessche P., Kirschner D. and Yakubu A.-A. (2002), Mathematical approaches for emerging and reemerging infectious diseases: models, methods, and theory, Vol. 126, Springer Science & Business Me- dia. Castillo-Chavez C. and Song B. (2004), ‘Dy- namical models of tuberculosis and their applications’, Mathematical Biosciences & Engineering 1(2), 361. COVID-19 pandemic (2022), ‘Covid-19 pan- demic in ethiopia — Wikipedia, the free encyclopedia’. Diekmann O., Heesterbeek J. and Roberts M. G. (2010), ‘The construction of next- generation matrices for compartmental epidemic models’, Journal of the Royal Society Interface 7(47), 873–885. Fleming W. H. and Rishel R. W. (2012), De- terministic and stochastic optimal control, Vol. 1, Springer Science & Business Me- dia. Gurmu E. D., Batu, G. B. and Wameko, M. S. (2020), ‘Mathematical model of novel covid-19 and its transmission dy- namics’, International Journal of Mathe- matical Modelling & Computations 10(2 (SPRING)), 141–159. Kifle Z. S. and Obsu, L. L. (2022), ‘Math- ematical modeling for covid-19 transmis- sion dynamics: A case study in ethiopia’, Results in Physics p. 105191. Lemecha Obsu L. and Feyissa Balcha S. (2020), ‘Optimal control strategies for the transmission risk of covid-19’, Journal of biological dynamics 14(1), 590–607. Lenhart S. and Workman J. T. (2007), Op- timal control applied to biological models, CRC press. Martcheva M. (2015), An introduction to mathematical epidemiology, Vol. 61, Springer. Pontryagin L. S. (2018), Mathematical theory of optimal processes, Routledge. Sohrabi C., Alsafi Z., O’Neill N., Khan,M., Kerwan A., Al-Jabir A., Iosifidis C. and Agha R. (2020), ‘World health organiza- tion declares global emergency: A review of the 2019 novel coronavirus (covid-19)’, International journal of surgery 76, 71– 76. Toquero C. M. (2020), ‘Challenges and op- portunities for higher education amid the covid-19 pandemic: The philippine con- text.’, Pedagogical Research 5(4). Valcher M. E. (2002), ‘Positive systems in the behavioral approach: main issues and recent results’, Electronic proceedings of MTNS 2002. Yang C. and Wang J. (2020), ‘A mathe- matical model for the novel coronavirus epidemic in wuhan, china’, Mathemat- ical biosciences and engineering: MBE 17(3), 2708. East Afr. J. Biophys. Comput. Sci. (2024), Vol. 5, No. 2, 58-86 86 EAJBCS__Volume_5_Issue 2_ Article 5_ Tinaw and Kiros_ A mathematical model for the transmission dynamics of COVID 19 pandemic.pdf (p.64-92) INTRODUCTION MATHEMATICAL MODEL FORMULATION The modified mathematical model Model formulation QUALITATIVE ANALYSIS OF THE MODIFIED MODEL Well-posedness Steady state Basic reproduction number Local stability of disease free equilibrium Global stability of disease free equilibrium Endemic equilibrium point Local stability of endemic equilibrium Bifurcation analysis Sensitivity analysis NUMERICAL SIMULATIONS AND DISCUSSION EXTENSION OF THE MODIFIED MODEL INTO AN OPTIMAL CONTROL Optimal protection and hospitalization using modified model Existence of an optimal control The Hamiltonian and optimality system Numerical simulations of optimal control problem Optimal control comparisons and strategies CONCLUSION