






























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 − Ĝ(Z1, Z2), Ĝ(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,

Ĝ(Z1, Z2) = AZ2 −G(Z1, Z2).

After some simplification, we obtain

Ĝ(Z1, Z2) =

 βσ2E(S
0 − S) + βσ1I(S

0 − S) + βσ3H(S0 − S)
0
0

 ,
Ĝ(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


