EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
2025, Vol. 18, Issue 4, Article Number 6596
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

Stability and Optimal Control Analysis of an SEIQR
Epidemic Model with Saturated Incidence Rate

Noshi Gul1, Ismail Shah2, Saeed Ahmad1,∗, Ihsan Ullah1, Manuel De la Sen3

1 Department of Mathematics, University of Malakand, Chakdara, Dir(L), 18800, Pakistan
2 Department of Mathematics, University of Nottingham Ningbo China, 199 Taikang East
Road, Ningbo 315100, China
3 Department of Electricity and Electronics, Institute of Research and Development of
Process, University of the Basque Country, Campus of Leioa (Bizkaia), Leioa 48940, Spain

Abstract. In this article, we proposed a new mathematical model to investigate the dynamics
of the infectious disease, control, and general disease transmission. The model exhibits two dis-
tinct non-trivial equilibrium states. As a fundamental prerequisite for stability analysis, we first
derive the epidemiological threshold parameter R0 through next-generation matrix methodology.
According to our investigation, R0 plays an essential role in describing the model’s dynamics.
We demonstrate that in the case when R0 takes values less or greater than unity, the endemic
(disease-free) condition is asymptotically stable both locally and globally. To try to stop the gen-
eral disease from spreading throughout a community, we add control parameters, create a control
model, and suggest control techniques. The maximum principle of Pontryagin is used to derive
the optimality system. Finally, the numerical simulations are performed using the fourth-order
Runge-Kutta technique to validate and confirm our analytical conclusions. Phase portrait analysis
further illustrates the convergence of system trajectories toward disease-free or endemic equilibria
under different control scenarios, reinforcing the stability criteria derived for R0.
2020 Mathematics Subject Classifications: 92D30, 93C15, 93D20, 34D23
Key Words and Phrases: Mathematical model, epidemic disease, asymptotical stability, optimal
control theory, phase portrait, numerical simulation

1. Introduction

Infectious diseases are one of the main causes of mortality around the globe. Infectious
diseases have existed for as long as humans have been on earth. Due to the emergence
and reemergence of some catastrophic diseases in recent decades, infectious diseases have
attracted the attention of many researchers [1–5] The main challenges regarding infectious

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6596

Email addresses: noshiguluom@gmail.com (N. Gul),
ismail8leu@gmail.com (I. Shah), saeedahmad@uom.edu.pk (S. Ahmad), ihsansf3@gmail.com (I. Ullah),
manuel.delasen@ehu.eus (M. De la Sen)

https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0)



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 2 of 23

diseases include studying the nature of their spread by considering possible factors causing
the diseases and predicting their future states. One way is to express these relations in
mathematical language. To look into it the mechanisms of infectious disease transfer
and forecast its future spread, mathematical models are crucial [6]. The estimation of
parameters, the development and testing of hypotheses that aid in disease prediction and
disease control are some of the fundamental characteristics of mathematical models. The
health department can use the mathematical model’s valuable guidelines to take action
towards the control and eradication of diseases, which is one of its most significant features.
Such models are helpful for identifying how people are simultaneously infected with a
disease and what strategy should be used to treat the condition brilliantly. Mathematical
modelling is one of the effective methods for predicting the dynamics of infectious diseases
in the realm of applied sciences. Numerous publications have been read to research the
dynamics and control of numerous infectious diseases [7–9] using mathematical modelling
of infectious diseases, which has a rich literature [10–15]. In order to research how to
control the transmission of viruses, some modified epidemic models are developed due
to the strong resemblance between software viruses and viruses that are living. The
results of the investigation have shown how crucial the incidence rate is to understanding
the nature of epidemic models [16]. Mathematical modeling uses various incident rates
to study the transmission dynamics of infectious disease [17]. Let h denote the disease
transmission rate, S and I are respectively used for susceptible and infected individuals.
The ratio hSI denotes an incidence rate, known as the bilinear incidence rate. This type
of incidence rate is used in various epidemic models [18–21]. The saturated incidence rate
hSI
1+αI was first proposed by Capasso and Serlo [22]. Whereas hI

1+αI tends to a saturation
level when I is big hI measures the force of infection when the disease is entering a fully
susceptible population, and 1

1+αI measures the inhibition effect from behavioural changes
in susceptible individuals as their numbers rise or from the overpopulation effect of the
infected people. Due to the inclusion of the behavioral changes and crowding effects of the
infected peoples and the prevention of the contact rate’s unboundedness by the selection
of acceptable parameters, this incidence rate is more rational than the bilinear incidence
rate and has been utilized in numerous pandemic simulations.

2. Model Formulation

In this section, we construct a general disease model and partition the total population
N(t) into five parts: S specially the susceptible individuals, E the Exposed people, I
the infectious density, Q the quarantine density, and R the people which recovered from
disease. The saturated incidence rate hSI

1+αI is adopted to enhance biological realism. Unlike
the bilinear form hSI, it incorporates saturation effects: as the infected population I grows,
the term 1

1+αI models reduced transmission due to behavioral changes or limited contact
capacity. This prevents unrealistic unbounded growth in infection force, offering a more
accurate representation of real-world epidemic dynamics. The suggested general disease



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 3 of 23

model is represented by the set of autonomous differential equations below:

dS

dt
= g − hSI

1 + αI
− (m+ b1)S,

dE

dt
=

hSI

1 + αI
− (m+ c+ b2)E,

dI

dt
= cE − (m+ b3 + q + µ) I,

dQ

dt
= b1S + b2E + b3I − (m+ p)Q,

dR

dt
= qI + pQ−mR.

(1)

The system indicated above is subject to the following preconditions. S(0) > 0, E(0) ≥ 0,
I(0) ≥ 0, Q(0) ≥ 0, R(0) ≥ 0.
N(t) = S + E + I +Q+R displays the total size of the community.
With the use of the derivative and model (3.0.1), we can

dN

dt
= g −mS −mE −mI −mQ−mR− µI = g −m(S + E + I +Q+R)− µI.

However, S + E + I +Q+R = N, the final equation has the following form:

dN

dt
= g −mN − µI.

This may be stated even more simply as:

dN

dt
≤ g −mN. (2)

The final inequality can be solved as follows by integrating with respect to t and using
the provided initial conditions:

N(t) =
g

m
+ Ce−mt.

Simplifying further and using the initial conditions N(0) = N0, we obtain

N(t) ≤ g

m
+

[
N(0)− g

m

]
e−mt.

On re-arranging we have,

N(t) ≤ g

m
(1− e−mt) +N(0)e−mt

It is confirmed by initial conditions that N(0) ≥ 0. It may also be noted from the last
inequality that the feasible area of the system is reached when the total population N(t)
remains positive and bounded.



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 4 of 23

3. Positivity of Solution

This proves that every phase pathway that begins in the positive area R5 ultimately
occurs in the phase space reaches and remains in the possibility area triangle (the model
falls into the biologically possible area). This can be accomplished by demonstrating that
the triangle is a positively invariant set and the systems global attractor. The positivity
of the proof is determined by the following theorem.

Theorem 1. The system (1) novel general disease model has a positive solution for all
beginning values that are supplied.

Proof. The first equation in the model, may be extracted as follows.
dS

dt
= g − hSI

1 + αI
− (m+ b1)S.

This can be written as:
dS

dt
≥ (λ(t) + (m+ b1))S ≥ g,

By using the given initial conditions, this first order differential equation’s solution is
provided as:

S(t) ≥ e−(m+b1)t−
∫ t
0 λ(t)dt

(
S(0) + g

∫ t

0
[e(m+b1)t1+

∫ t
0 λ(t)dt]

)
.

Due to the fact that under integration all the constants in the suggested model are positive
which possesses the nonnegativity in both the invariant and results:

△ =

{
S(t), E(t), I(t), Q(t), R(t) ∈ R5 : S > ‘0(E, I,Q,R ≥ 0), N(t) ≤ g

m

}
.

4. Steady States

To determine how the suggested model behaves qualitatively, there are two different
kinds of stable states: endemic equilibriums and equilibriums absence of illness. The
endemic equilibrium point is represented by E1, while the equilibrium point devoid of
illness is represented by E0.

4.1. Disease Free Equilibrium point, E0

The disease-free stationary point designates the location at which the illness has been
totally eradicated from the local population. For the suggested model’s no-disease equi-
librium to be found, given that the illness is not found in the population, every equation
has a zero on its right hand, resulting in E = I = 0. The no-sickness stationary point is
achieved by computation and is denoted by E0, as seen below:

E0 = (S0, E0, I0, Q0, R0) =

(
g

m+ b1
, 0, 0,

b1g

(m+ p)(m+ b1)
,

b1gp

m(m+ p)

)
.



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 5 of 23

4.2. Endemic Equilibria
According to the endemic equilibrium, this occurs when the illness affects a community

for a lengthy period of time. We shall put the right side of the proposed model equal to
zero in order to locate the endemic equilibria. The endemic equilibria for the system
(3.0.1) at E1 = (S∗, E∗, I∗, Q∗, R∗) are

S∗ =
(m+ c+ b2)(m+ b3 + q + µ)(1 + αI∗)

hc
,

E∗ =
(m+ b3 + q + µ)I∗

c
,

I∗ =
hgc− (m+ b1)(m+ b2 + c)(m+ b3 + q + µ)

(m+ b2 + c)(m+ b3 + q + µ)(h+ α(m+ b1))
,

Q∗ =
b1S

∗ + b2E
∗ + b3I

∗

m+ p
,

R∗ =
qI∗ + pQ∗

m
.

4.3. Basic Threshold Number
The basic reproduction number R0 is defined as the average number of secondary in-

fections produced by one infected individual in a completely susceptible population. This
time, R0 is calculated as typically there are how many secondary cases of general ill-
ness peeople that an individual with a general disease who was not receiving treatment
throughout his or her general disease caused in a population of potential general disease
patients.
The well-known next generation matrix may be used to establish the fundamental repro-
duction number. For this, we separate patients who are infected and those who are not.
According to model, the infected classes are E and I. Consider only those classes which
infection of the disease posses, the model becomes

J1 = (E, I),

dJ1
dt

=
d

dt
(E, I).

dE

dt
=

hSI

1 + αI
− (m+ c+ b2)E,

dI

dt
= cE − (m+ b3 + q + µ)I. (3)

Using the system’s making the matrices F and V which show the susceptible people and
infected people.

dx

dt
= F − V ,



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 6 of 23

Where,

F =

[
hSI
1+αI

0

]
, V =

[
(m+ c+ b2)E

−cE + (m+ b3 + q + µ)I

]
,

Then by finding the maximum spectram of the system F ∗V −1∗

R0 =
hcg

(m+ b1)(m+ c+ b2)(m+ b3 + q + µ)
.

which is the basic reproductive value R0.
The basic reproduction number R0 acts as a key epidemiological threshold governing

system dynamics. If R0 < 1, the disease-free equilibrium is both locally and globally
asymptotically stable, leading to disease eradication. If R0 > 1, the disease-free equilib-
rium becomes unstable, and the endemic equilibrium is asymptotically stable, reflecting
disease persistence. Thus, R0 directly determines the long-term behavior of the infection
within the population.



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 7 of 23

0

1

0.2

0.4

1

R
0

0.6

0.8

b
2

0.8

0.5 0.6

b
1

1

0.4
0.2

0 0

0

1

2

1

R
0

4

0.8

µ

0.5 0.6

b
2

6

0.4
0.2

0 0

1

1

2

1

R
0

0.8

3

µ

0.5 0.6

b
3

0.4
0.2

0 0

0

1

2

1

R
0

4

0.8

µ

0.5

6

0.6

c

0.4
0.2

0 0

0

1

20

1

R
0

40

0.8

m

0.5 0.6

h

60

0.4
0.2

0 0

0

1

1

1

R
0

2

0.8

µ

0.5 0.6

b
1

3

0.4
0.2

0 0

Figure 1: Behaviour of the basic reproduction number in three dimensions in relation to
different model parameters.



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 8 of 23

5. Local Stability Analysis

In this section, we examine the model’s (1) endemic equilibrium point (E∗) and without
illness equilibrium point’s (E0) local and global stability analyses.

5.1. Local Stability Analysis of General Disease Model
Theorem 2. If R0 < 1, then the system’s disease-free equilibrium point E0 is locally
asymptotically stable.

Proof. The Jacobin matrix at disease free equilibrium is;

J(E0) =


−(m+ b1) 0 − hg

m+b1
0 0

0 −(m+ c+ b2)
hg

m+b1
0 0

0 c −(m+ b3 + q + µ) 0 0
b1 b2 b3 −(m+ p) 0
0 0 q p −m

 .

The three eigenvalues are λ1 = −m, λ2 = −(m+ p) and λ3 = −(m+ b1). To find the next
eigenvalues, we use the 2× 2 matrix and after algebraic manipulation we reached to

λ2 + λ(2m+ c+ b2 + b3 + q + µ) + (m+ c+ b2)(m+ b3 + q + µ)− hgc

m+ b1
= 0. (4)

When R0 > 1, then the above equation (m+c+b2)(m+b3+q+µ)− hgc
m+b1

< 0, This implies
that Eq. (4) has a positive and negative root. Therefore, the disease-free equilibrium E0

is unstable.

5.2. Global Asymptotic Stability of the Disease Free Equilibrium
In order to demonstrate the globally asymptotically stable nature of the disease-free

equilibrium E0, we examine the Lyapunov function for R0 ≤ 1,.

H(E, I) = B1E +B2I

The derivation of H(E, I) with regrd to t gives,

dH

dt
= B1

dE

dt
+B2

dI

dt
,

dI

dt
= B1

( hSI

1 + αI
− (m+ c+ b2)E

)
+B2

(
cE − (m+ b3 + q + µ)I

)
,

= B1
hSI

1 + αI
−B1(m+ c+ b2)E +B2cE −B2(m+ b3 + q + µ)I,

(5)



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 9 of 23

calculating the above equation which gives B1 = c, and B2 = (m+ c+ b2). so,

dH

dt
=

[
chS

1 + αI
− (m+ c+ b2)(m+ b3 + q + µ)

]
≤

[
chS

1 + αI
− (m+ c+ b2)(m+ b3 + q + µ)

]
≤ chS0

R0

[
R0S

S0
− 1

]
I ≤ 0.

(6)

Furthermore dH(E,I)
dt = 0 if and only if I = 0. In light of this, the greatest compact

invariant set in
{
(S,E, I,Q,R)|H(E, I) = 0

}
. When R0 ≤ 1, is the singleton E0. Lasalle’s

invariance principle implies that E0 is locally asymptotically stable.

5.3. Local Asymptotic Stability of Endemic Equilibrium
Theorem 3. If 1 < R0 then (1) at endemic equilibrium, is locally asymptotically stable.
The situation of endemicity state is unstable whenever 1 > R0.

Proof. Jacobian matrix J at E∗ shows that

J (E∗) =


− hI∗

1+αI∗ − (m+ b1) 0 hS∗

(1+αI∗)2 0 0
hI∗

1+αI∗ −(m+ c+ b2)
hS∗

(1+αI∗)2 0 0

0 c −(m+ b3 + q + µ) 0 0
b1 b2 b3 −(m+ p) 0
0 0 q p −m

 .

Eigenvalues of J(E∗) are negative numbers that are λ1 = −m, λ2 = −(m + p) then we
only need to consider the roots of λ3 +B1λ

2 +B2λ+B3 = 0, where

B1 =
hI∗

1 + αI∗
+ 3m+ 2b1 + b3 + c+ q > 0,

B2 =

(
hI∗

1 + αI∗
+m+ b1

)
(2m+ b2 + b3 + c+ q + µ) + (m+ c+ b1) (m+ b3 + q + µ)− chS∗

(1 + αI∗)2

>

(
hI∗

1 + αI∗
+m+ b1

)
(2m+ b2 + b3 + c+ q + µ) > 0,

B3 =

(
hI∗

1 + αI∗
+m+ b1

)(
(m+ c+ b1) (m+ b3 + q + µ)− chS∗

(1 + αI∗)2

)
+

ch2S∗I∗

(1 + αI∗)3
> 0.

Based on the following relation chS∗

1+αI∗ = (m+ c+ b1) (m+ b3 + q + µ) .
By a direct calculation, we have that B1B2−B3 > 0. Following Routh-harwitz criteria one
can easily confirm endemic equilibrium E∗ is locally asymptotically stable. This completes
the proof.



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 10 of 23

6. Global Stability Analysis

Theorem 4. If R0 > 1, the endemic equilibrium E∗ is globally asymptotically stable.

Proof. Taking the subsystems in this case,

dS

dt
= g − hSI

1 + αI
− (m+ b1)S,

dE

dt
=

hSI

1 + αI
− (m+ c+ b2)E,

dI

dt
= cE − (m+ b3 + q + µ) I. (7)

The Jacobian matrix of the system is given as J.

J =

 − hI
1+αI − (m+ b1) 0 − hs

(1+α)2
hI

1+αI − (m+ c+ b2)
hs

(1+α)2

0 c −(m+ b3 + q + µ)

 ,

where J |2| stands for the second additive matrix of J. This matrix form exists in

J |2| =

 − hI
1+αI − k hS

(1+αI)2
hS

(1+α)2

c − hI
1+αI − l 0

0 hI
1+αI −m

 .

Entries of the last matrix are defined in the following way,
k = 2m+ c+ b1 + b2,
l = 2m+ b1 + b3 + q + µ,
m = 2m+ c+ b2 + b3 + q + µ, Let us consider the function

X(x) = X (S,E, I) = diag

{
1,

E

I
,
E

I

}
Then,

X−1 = diag

{
1,

I

E
,
I

E

}
, Xf = diag

{
0,

Ė

I
− Eİ

I2
,
Ė

I
− Eİ

I2

}
.

Direct calculation demonstrates that which can be seen further as

XfX
−1 = diag

{
0,

Ė

E
− İ

I
,
Ė

E
− İ

I

}
.

XJ |2|X−1 =

 x11
hSI

E(1+αI)2
hSI

E(1+αI)2
Ec
I x22 0

0 hI
1+αI x33





N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 11 of 23

Where
x11 =

−hI
1+αI − (2m+ c+ b1 + b2),

x22 =
−hI
1+αI − (2m+ b1 + b3 + q + µ),

x33 = −(2m+ c+ b2 + b3 + q + µ).
Thus, we are able to write

Q = XfX
−1 +XJ |2|X−1.

It is simple to check if the matrix B has the specified form after values have been inserted.

Q =

[
Q11 Q12

Q21 Q22

]
.

Take note that the entries in matrix Q are computed as

Q11 = x11,

Q12 =
hSI

E(1 + αI)2
,

hSI

E(1 + αI)2
,

Q21 =

[
cE
I
0

]
,

Q22 =

[
x22 +

IĖ−Eİ
IE 0

hI
1+αI x33 +

IĖ−Eİ
IE

]
.

Let (u, v, w) be a vector in R3. Its norm ∥·∥ is defined as

∥(u, v, w)∥ = max{|u|, |v|+ |w|}.

Let µ(Q) be the Lozinskiĭ measure with respect to this norm. Then

µ(Q) ≤ sup{g1, g2},

where

g1 = µ1(Q11) + ∥Q12∥,
g2 = ∥Q21∥+ µ1(Q22).

Here ∥Q12∥ and ∥Q21∥ are matrix norms with respect to the ℓ1 vector norm, and µ1 denotes
the Lozinskiĭ measure with respect to this ℓ1 norm. Then

µ1(Q11) = − hI

1 + αI
− (2m+ c+ b1 + b2),

∥Q21∥ =
cE

I
,

∥Q12∥ = max

{
hSI

E(1 + αI)2
,

hSI

E(1 + αI)2

}
=

hSI

E(1 + αI)2
,



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 12 of 23

and

µ1(Q22) = max

{
Ė

E
− İ

I
− hI

1 + αI
−N,

hI

1 + αI
,

Ė

E
− İ

I
−K

}
=

Ė

E
− İ

I
−min(N,K).

Therefore, we have

g1 = − hI

1 + αI
−M +

hSI

E(1 + αI)2
,

g2 =
cE

I
+

Ė

E
− İ

I
−min(N,K).

From the proposed model (1.1), we get

Ė

E
=

hSI

E(1 + αI)
− (m+ c+ b2),

İ

I
=

cE

I
− (m+ b3 + q + µ).

Then we have

g1 = − hI

1 + αI
− (2m+ c+ b1 + b2) +

hSI

E(1 + αI)2
+

[
hSI

E(1 + αI)2
− hSI

E(1 + αI)

]
≤ − hI

1 + αI
− (m+ b1) +

Ė

E
,

g2 =
cE

I
+

Ė

E
− İ

I
−min(N,K)

≤ Ė

E
−m.

Since
cE

I
=

İ

I
− (m+ b3 + q).

Furthermore, we obtain:

µ(Q) ≤ sup{g1, g2} ≤ sup

{
− hI

1 + αI
− (m+ b1) +

Ė

E
,
Ė

E
−m

}
≤ Ė

E
−m,

then
1

t

∫ t

0
µ(Q) ds ≤ 1

t

∫ t

0

(
Ė

E
−m

)
ds =

1

t
ln

E(t)

E(0)
−m,

which implies q ≤ m
2 < 0. We know that positive equilibrium (S∗, E∗, I∗) is globally

asymptotically stable.
In this section, it is described how the general disease model may be used to optimise

control concepts.



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 13 of 23

7. Formulation of Optimal Control Problem

In order to design our control approach in this part, we take into consideration the
general illness pandemic model. Through the use of the three control variables u1(t),
which denote the effort that reduces the contact between the susceptibles and infectious
individuals; u2(t), denote the rate at which infectious individuals are treated and u3(t),
shows the vaccination coverage. Thus, by including the above-mentioned control variables,
the general disease control model is created, and it looks like this:

dS

dt
= g − (1− u1)

hSI

1 + αI
− (m+ b1)S − u3S,

dE

dt
= (1− u1)

hSI

1 + αI
− (m+ c+ b2)E,

dI

dt
= cE − (m+ b3 + µ) I − (1 + u2) qI,

dQ

dt
= b1S + b2E + b3I − (m+ p)Q,

dR

dt
= (1 + u2)qI + pQ−mR+ u3S.

(8)

The system (8) is subject to the nonnegative initial conditions [23, 24]

S(0) > 0, E ≥ 0, I ≥ 0, Q ≥ 0 and R(0) ≥ 0.

We define the goal functional function, which aims to decrease the spread of the general
disease infection, can be described as;

J(u1(t), u2(t), u3(t)) =

∫ Tf

0
A1E +A2I −A3R+

1

2

[
A4u

2
1(t) +A5u

2
2(t) +A6u

2
3(t)

]
.(9)

The values A1, A2 and A3 are weight constants. The objective functional serves the pur-
pose of minimising the number of affected individuals. In this way, we determine an
optimal control triplet:

J(u∗1(t), u
∗
2(t), u

∗
3(t)) = min{J(u1(t), u2(t), u3(t)) ∈ U}. (10)

According to the problem being considered, the control set is provided by
U = {u1, u2, u3) : [0, Tf ] → [0, 1], (u1, u2, u3)is a Labesgue measurable}.
Next, we demonstrate that the optimum control problem exists.

7.1. Existence of the Optimal Control Problem
Here, we provide proof that shows the control problem exists. For the optimum control

issue, we provide the following definition of Hamiltonian H:

H = L(E(t), I(t), R(t), u1(t), u2(t), u3(t)) + λ1
dS

dt
+ λ2

dE

dt
+ λ3

dI

dt
+ λ4

dQ

dt
+ λ5

dR

dt
.(11)



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 14 of 23

Theorem 5. For control issue (4.1), there is

u∗(t) = (u∗1(t), u
∗
2(t), u

∗
3(t)) ∈ U,

in such away that

min(u1(t),u2(t),u3(t)∈U)J(u1(t), u2(t), u3(t)) = J(u∗1, u
∗
2, u

∗
3(t)).

Proof. For the purpose of validating the optimal control rate, we used a variety of
methodologies. Thus, each and every control and state variable has a positive value. Due to
this, the issue is being reduced, So, u1(t), u2(t) and u3(t) elaborates the required convexity
of the goal functional and is satisfied. The collection of control variables u1, u2, u3 ∈ U
is by definition convex and closed. There is a defined optimum system, Furthermore, it
supplies the confidence regarding the solidity required for the verification of the optimal
control system. Moreover, an integral in the practical purpose A1E+A2I−A3R+ 1

2(A4u
2
1+

A5u
2
2 +A6u

2
3) which confirms the evidence, is convex on the control set U.

8. Optimality System

Theorem 6. Given optimal controls u∗1(t), u
∗
2(t), u

∗
3(t) and solutions S∗, E∗, I∗, Q∗, R∗ of

the corresponding state system (8) there is adjoint variables λn(t), n = 1, ..., 5,

dλ1

dt
= (λ1 − λ2)(1− u1)

hI

1 + αI
+ (λ1 − λ4)b1 + (λ1 − λ5)u3 +mλ1,

dλ2

dt
= −A1 + (λ2 − λ3)c+ (λ2 − λ4)b2 +mλ2,

dλ3

dt
= −A2 + (λ1 − λ2)(1− u1)

hS

(1 + αI)2
+ (λ3 − λ4)b3 + (λ3 − λ5)(1 + u2)q + (m+ µ)λ3

dλ4

dt
= (λ4 − λ5)p+mλ4

dλ5

dt
= −A3 +mλ5,

with transversality conditions λn(Tf ) = 0, n = 1, ..., 5.

Proof. When we consider the values as S(t) = S∗, E(t) = E∗, I(t) = I∗, Q(t) =
Q∗, and R(t) = R∗ and differentiate the Hamiltonian with respect to state variables
S(t), E(t), I(t), Q(t) and R(t), respectively, we get the adjoint system (12) with transver-
sality conditions λn(t) = 0, n = 1, 2, ..., 5. In coming part, we develop the principles of
optimal control variables as follows: while using optimal conditions for solution ∂H

∂u1
=

0, ∂H
∂u2

= 0 and ∂H
∂u3

= 0 on the inside of the control setup.
Lastly, we employ the control area characteristic while writing.

u1(t) =
(λ2 − λ1)

A4

hSI

1 + αI
,



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 15 of 23

u2(t) =
(λ3 − λ5)qI

A5
,

u3(t) =
(λ1 − λ5)S

A6
. (12)

As one can see, the equation (12) for U∗ = (u∗1, u
∗
2, u

∗
3) is made reference to as a represen-

tation of optimal controls. State variables and optimum control are obtained by solving
the optimality system. When values are entered for u1, u2 and u3 in the system, we obtain
the underlying syste,:

dS∗

dt
= g − hS∗I∗

1 + αI∗
− (m+ b1)S

∗ − (
(λ1 − λ5)S

A6
)S∗,

dE∗

dt
= (1− (

(λ2 − λ1)

A4

hSI

1 + αI
))

hS∗I∗

1 + αI∗
− (m+ c+ b2)E

∗,

dI∗

dt
= cE∗ −

(
m+ b3 + (1 + (

(λ3 − λ5)qI

A5
))q + µ

)
I∗,

dQ∗

dt
= b1S

∗ + b2E
∗ + b3I

∗ − (m+ p)Q∗,

dR∗

dt
= (1 + (

(λ3 − λ5)qI

A5
))qI∗ + pQ∗ −mR∗ + (

(λ1 − λ5)S

A6
)S∗.

(13)

The next section includes several numerical simulations. To verify the analytical find-
ings from the previous section. The numerical simulations are also briefly discussed.

9. Numerical Results

9.1. Numerical Simulations for the Problem (Without Control Variables)
The numerical simulations used to validate the results of analysis examined in the ear-

lier chapters are the focus of this portion of the thesis. We employ the widely established
Runge-Kutta method of fourth order (RK4). In the beginning, we run numerical simu-
lations for a case in which an area is free of disease, i.e the so-called condition devoid of
infections. The threshold quantity’s value is assumed to be smaller than unity by selecting
parameter values in the model under consideration.



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 16 of 23

0 20 40 60 80 100

Time 

20

40

60

80

100

120

140

160

180

S
u

sc
ep

ti
b

le
 C

la
ss

S(0)=180

S(0)=160

S(0)=140

S(0)=120

0 20 40 60 80 100

Time 

0

20

40

60

80

100

120

140

160

180

E
xp

o
se

d
 C

la
ss

E(0)=180

E(0)=160

E(0)=140

E(0)=120

0 50 100 150

Time 

0

50

100

150

200

250

300

In
fe

ct
ed

 C
la

ss

I(0)=180

I(0)=160

I(0)=140

I(0)=120

0 50 100 150

Time 

0

50

100

150

200

250

300

Q
u

a
ra

n
ti

n
e
 C

la
ss

Q(0)=180

Q(0)=160

Q(0)=140

Q(0)=120

0 20 40 60 80 100

Time 

100

150

200

250

300

350

400

450

500

R
e
c
o

v
e
r
e
d

 C
la

s
s

R(0)=180

R(0)=160

R(0)=140

R(0)=120

Figure 2: Dynamical behaviour of the trajectories of the total classes of the model for the
no-infection state.

In Fig. (2) (Susceptible population graph) we plot the trajectories of the susceptible
compartment of the model under consideration for the infection-free condition versus time
t, for varied initial sizes of the compartment. From this figure, we observe that the
susceptible population decreases and gradually tend to the no-disease state S0. We also
observe that this equilibrium is nonzero. It means that there will always be chance to
people to catch the disease, which is biologically relevant. (Exposed population graph)
demonstrates how the exposed class of infection behave dynamically. The class is rapidly
decreasing at the no illness state. As time increases, the class becomes closer to the zero
state, or the equilibrium free of sickness. In Fig. (Infected population graph) we depict
the infected class. At the disease-free condition, there is a quickly increase for first few
weeks in the class. As time increases, the class becomes closer to the zero condition,



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 17 of 23

or the equilibrium free from sickness. In other words, the infection population reaches
its convergence at the zero condition as time goes on. This agrees with the value of I0

already determined for the no-infection state. In Fig. (Quarantine population graph)
the quarantine class’s solution curve is displayed. The class is growing over time and
eventually converges to an infection-free state, demonstrating that this class’s equilibrium
is non-zero. In Fig. (Recovered population graph) different beginning conditions are used
to run the simulations. As time goes on, we see a sharp rise in the recovered population
and a tendency of the class towards zero, the no-infection equilibrium condition. Now,
we consider the case of persistence of the disease in a community. For this we take the
values: g = 10, h = 0.01, m = 0.0124, b1 = 0.05, c = 0.04, b2 = 0.0246, b3 = 0.015, q =
0.2, µ = 0.2, p = 0.4, α = 0.2 and perform numerical simulations.

0 20 40 60 80 100

Time 

90

100

110

120

130

140

150

160

170

180

S
u

sc
ep

ti
b

le
 C

la
ss

S(0)=180

S(0)=160

S(0)=140

S(0)=120

0 10 20 30 40 50

Time 

0

20

40

60

80

100

120

140

160

180

E
xp

o
se

d
 C

la
ss

E(0)=180

E(0)=160

E(0)=140

E(0)=120

0 10 20 30 40 50

Time 

0

50

100

150

200

250

300

In
fe

ct
ed

 C
la

ss

I(0)=180

I(0)=160

I(0)=140

I(0)=120

0 10 20 30 40 50

Time 

0

20

40

60

80

100

120

140

160

180

Q
u

a
ra

n
ti

n
e
 C

la
ss

Q(0)=180

Q(0)=160

Q(0)=140

Q(0)=120

0 100 200 300 400 500

Time 

0

500

1000

1500

2000

2500

3000

R
e
c
o

v
e
r
e
d

 C
la

s
s

R(0)=180

R(0)=160

R(0)=140

R(0)=120

Figure 3: Dynamical behaviour of the trajectories of the total classes of the model for the
infection state.



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 18 of 23

Fig. (3) displays the infected compartment in the event that the disease continues to
spread throughout a community for different beginning compartment sizes. The compart-
ment is initially shown to be rising; as time goes on, the compartment falls and converges
to the illness endamic state. This validates our analytical findings. While in the last figure
(recovered population graph) show the picture of the recovered population versus time for
various initial conditions. We see that throughout the first several weeks, this class grows
and then progressively approaches the endemic stability threshold.

9.2. Numerical Simulations for the Optimal Control Problem
We numerically solve the optimal control model (8) using the Runge-Kutta method of

order 4. The parametric values are taken as: A1 = 2, , A2 = 3, A3 = 4, A4 = 5, A5 =
A6 = 1, g = 10, h = 0.00245, m = 0.0124, b1 = 0.4, c = 0.4, b2 = 0.5, b3 = 0.45, q =
0.5, µ = 0.023, p = 0.4, α = 0.2, The Fig. (4) are the solution curves of the susceptible
class both in the absence and presence of the optimal control parameters. We observe from
the figure, that the susceptible population can be increased once one applies the control
stratifies. In other words, the number of susceptible individuals increase when the control
measures are adopted. The figure also show that the class S(t) decreases whenever these
measures are avoided.

0 10 20 30 40 50

Time

0

200

400

600

800

1000

S
u

sc
e
p

ib
le

 C
la

ss

With control

Without control

Figure 4: Solution curves of the trajectories of susceptible class of the model under con-
sideration with and without control parameters.

Plotted in the second panel of Fig. (5) are the trajectories of the exposed and quar-
antine classes both in the absence and presence of the control measures. It is observed
from both figures that the population can be decreased with the application of the control
measures, which is the aim of the optimal control problem.



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 19 of 23

0 50 100 150 200 250 300

Time

0

100

200

300

400

500

600

700

E
x
p

o
se

d
 C

la
ss

Control in  the exposed class

Without control

With control

0 10 20 30 40 50

Time

0

500

1000

1500

C
o

n
tr

o
l 

o
n

 t
h

e
 c

la
ss

 Q
(t

)

Control in the Infectious class

Without control

With control

Figure 5: Solution curves of the trajectories of various classes of the model under consid-
eration with and without control parameters.

In the same way, other panels of the figure Fig. (6) shows that the infected popu-
lation can be decreased once the control measures are adopted. Similarly the recovered
populations are enhanced with the control measures suggested in our control problem.

0 10 20 30 40 50

Time

0

100

200

300

400

500

In
fe

ct
ed

 C
la

ss
 I

(t
)

without control

with control

0 50 100 150 200 250 300

Time

0

500

1000

1500

2000

R
e
c
o

v
e
re

d
 C

la
ss

 R
(t

)

Without control

With control

Figure 6: Solution curves of the trajectories of various classes of the model under consid-
eration with and without control parameters.

10. Phase Portrait Study for Proposed Model

The phase portrait analysis provides a comprehensive visualization of the system’s
trajectories under both controlled and uncontrolled scenarios, highlighting the impact of
intervention strategies on disease dynamics. System of figure (7) illustrates the interplay
between different populations, demonstrating distinct behavioral patterns. In the absence
of control measures, trajectories diverge toward higher infection levels, aligning with the
endemic equilibrium. Conversely, with optimal controls u1, u2, and u3 applied, trajectories



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 20 of 23

converge toward regions of reduced infection density, reflecting the efficacy of interventions
in steering the system closer to the disease-free equilibrium. Similarly, System of (8) also
compares the temporal evolution of different compartments, revealing that control imple-
mentation suppresses E while elevating Q, consistent with enhanced quarantine efforts.
The shrinking of phase space near endemic states under control shows stable dynamics,
conforming the global stability that has been mathematically proven. The visualizations
above show how control parameters influence the attractors of the system, proving their
role in reducing the severity of outbreaks. The convergence patterns we witness confirm
the theoretical stability criteria of R0, showing that efficient strategies can readily break
transmission networks, as testified by fewer cycles and faster stabilization. This geometric
analysis is corroborated by numerical results, providing easy-to-interpret evidence of the
model’s ability to make predictions under intervention.

0 50 100 150

S

0

1

2

3

4

5

6

7

E

S vs E (With/Without Control)

Without Controlled

With Controlled

Equilibrium

0 50 100 150

S

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

I

S vs I (With/Without Control)

Without Controlled

With Controlled

Equilibrium

0 50 100 150

S

0

2

4

6

8

10

12

14

Q

S vs Q (With/Without Control)

Without Controlled

With Controlled

Equilibrium

0 50 100 150

S

0

20

40

60

80

100

120

140

R

S vs R (With/Without Control)

Without Controlled

With Controlled

Equilibrium

0 1 2 3 4 5 6 7

E

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

I

E vs I (With/Without Control)

Without Controlled

With Controlled

Equilibrium

0 1 2 3 4 5 6 7

E

0

2

4

6

8

10

12

14

Q

E vs Q (With/Without Control)

Without Controlled

With Controlled

Equilibrium

0 1 2 3 4 5 6 7

E

0

20

40

60

80

100

120

140

R

E vs R (With/Without Control)

Without Controlled

With Controlled

Equilibrium

0 0.2 0.4 0.6 0.8 1

I

0

2

4

6

8

10

12

14

Q

I vs Q (With/Without Control)

Without Controlled

With Controlled

Equilibrium

0 0.2 0.4 0.6 0.8 1

I

0

20

40

60

80

100

120

140

R

I vs R (With/Without Control)

Without Controlled

With Controlled

Equilibrium

Figure 7: Phase Photrate system



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 21 of 23

0 2 4 6 8 10 12 14

Q

0

20

40

60

80

100

120

140

R

Q vs R (With/Without Control)

Without Controlled

With Controlled

Equilibrium

0 1 2 3 4 5 6 7

E

0

50

100

150

S

E vs S (With/Without Control)

Without Controlled

With Controlled

Equilibrium

0 0.2 0.4 0.6 0.8 1

I

0

50

100

150

S

I vs S (With/Without Control)

Without Controlled

With Controlled

Equilibrium

0 2 4 6 8 10 12 14

Q

0

50

100

150

S

Q vs S (With/Without Control)

Without Controlled

With Controlled

Equilibrium

0 20 40 60 80 100 120 140

R

0

50

100

150

S

R vs S (With/Without Control)

Without Controlled

With Controlled

Equilibrium

0 20 40 60 80 100 120 140

R

0

1

2

3

4

5

6

7

E

R vs E (With/Without Control)

Without Controlled

With Controlled

Equilibrium

0 2 4 6 8 10 12 14

Q

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

I

Q vs I (With/Without Control)

Without Controlled

With Controlled

Equilibrium

0 20 40 60 80 100 120 140

R

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

I

R vs I (With/Without Control)

Without Controlled

With Controlled

Equilibrium

0 20 40 60 80 100 120 140

R

0

2

4

6

8

10

12

14

Q

R vs Q (With/Without Control)

Without Controlled

With Controlled

Equilibrium

Figure 8: Phase Photrate system

11. Conclusions

We examined the mathematical model of epidemics in this manuscript. We talked
about a few key characteristics, such as the investigated model’s solutions’ positivity
and boundedness within a biologically acceptable range. We calculated the fundamental
threshold number R0 in order to talk about the stability analysis of the potential non-
trivial equilibria of the suggested model. This was accomplished by utilising a well-known
the next-generation techniques for technology. We’ve set criteria for both the local and
global stability of the simulation under consideration’s no-infection and disease-endemic
states. Using the linearization technique, we conducted the studies of local stability. The
Castillo-Chavez approach and geometric method are used accordingly to explore the global
stability assessments of a disease-endemic condition and the no-infection state. Research
has shown that when the fundamental threshold number R0 takes values smaller or larger
than unity, the no-infection (endemic) condition is maintained on a global and local scale.
In order to stop a general illness infection from spreading, we include three control vari-
ables: u1(t), which denote the effort that reduces the contact between the susceptible and
infectious individuals; u2(t), denote the rate at which infectious individuals are treated



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 22 of 23

and u3(t), shows the vaccination coverage. We numerically solve the proposed model us-
ing the widely recognised fourth order Runge-Kutta (RK-4) approach. Our numerically
obtained results well support the analytical findings. The phase portrait analysis further
demonstrates how system trajectories evolve under different control strategies, visually
confirming the stabilization toward disease-free equilibrium when interventions are ap-
plied. The current work may be extended to fractional order differential equations and
various approaches can be used to study the transmission dynamics of a general infection.
Work on such problems is under investigation and will be reported in a near publication.

Acknowledgements

Manuel De la Sen is thankful to the Basque Government, Grants IT1555-22, for the
financial support.

References

[1] A. Omame, H. Rwezaura, M. L. Diagne, S. C. Inyama, and J. M. Tchuenche. Covid-
19 and dengue co-infection in brazil: optimal control and cost-effectiveness analysis.
European Physical Journal Plus, 136(10):1090, 2021.

[2] A. Abidemi, Z. M. Zainuddin, and N. A. B. Aziz. Impact of control interventions on
covid-19 population dynamics in malaysia: a mathematical study. European Physical
Journal Plus, 136:1–35, 2021.

[3] K. Wang, A. Fan, and A. Torres. Global properties of an improved hepatitis b virus
model. Nonlinear Analysis: Real World Applications, 11(4):3131–3138, 2010.

[4] X. Liu, M. U. Rahman, S. Ahmad, D. Baleanu, and Y. N. Anjum. A new fractional
infectious disease model under the non-singular mittag–leffler derivative. Waves in
Random and Complex Media, pages 1–27, 2022.

[5] M. G. M. Gomes, A. Margheri, G. F. Medley, and C. Rebelo. Dynamical behaviour of
epidemiological models with sub-optimal immunity and nonlinear incidence. Journal
of Mathematical Biology, 51:414–430, 2005.

[6] O. Zakary, M. Rachik, and I. Elmouki. On the analysis of a multiregions discrete sir
epidemic model: an optimal control approach. International Journal of Dynamics
and Control, 5:917–930, 2017.

[7] S. Ullah, M. F. Khan, S. A. A. Shah, M. Farooq, M. A. Khan, and M. B. Mamat.
Optimal control analysis of vector-host model with saturated treatment. European
Physical Journal Plus, 135(10):1–25, 2020.

[8] G. Zaman, Y. H. Kang, and I. H. Jung. Stability and optimal vaccination of an sir
epidemic model. BioSystems, 93:240–249, 2008.

[9] A. V. Kamyad, R. Akbari, and A. Heydari. Mathematical modeling of transmission
dynamics and optimal control of vaccination and treatment for hepatitis b virus.
Computational and Mathematical Methods in Medicine, pages 1–15, 2004.

[10] G. Zaman, Y. H. Kang, and I. H. Jung. Optimal treatment of an sir epidemic model
with time delay. BioSystems, 98:43–50, 2009.



N. Gul et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6596 23 of 23

[11] I. Shah, H. Alrabaiah, and B. Ozdemir. Using advanced analysis together with frac-
tional order derivative to investigate a smoking tobacco cancer model. Results in
Physics, 51:106700, 2023.

[12] S. Ahmad, N. Ahmad, and I. Shah. Stability and sensitivity analysis of cyberattack
propagation models in computer networks. European Journal of Pure and Applied
Mathematics, 18(3):6336, 2025.

[13] I. Shah, I. Ali, A. Ali, I. Ahmad, S. Islam, G. Rasool, S. Formanova, and M. Kallel.
Optimal control and sensitivity analysis of a mathematical model for mdr-tb transmis-
sion with advanced treatment strategies. European Physical Journal Plus, 140(6):1–
15, 2025.

[14] I. Ullah, S. Ahmad, M. Arfan, and M. De la Sen. Investigation of fractional order
dynamics of tuberculosis under caputo operator. Fractal and Fractional, 7(4):300,
2023.

[15] H. Inaba. Mathematical analysis of an age-structured sir epidemic model with vertical
transmission. Discrete and Continuous Dynamical Systems - Series B, 6:69, 2006.

[16] X. Liu and L. Yang. Stability analysis of an seiqv epidemic model with saturated
incidence rate. Nonlinear Analysis: Real World Applications, 13(6):2671–2679, 2012.

[17] L. Esteva and M. Matias. A model for vector transmitted diseases with saturation
incidence. Journal of Biological Systems, 9(4):235–245, 2001.

[18] S. R. Chawla, S. Ahmad, and A. Khan. Sensitivity and bifurcation analysis of pine
wilt disease with harmonic mean type incidence rate. Physica Scripta, 97(5):055006,
2022.

[19] A. Korobeinikov and P. K. Maini. Non-linear incidence and stability of infectious
disease models. Mathematical Medicine and Biology, 22(2):113–128, 2005.

[20] M. D. Samsuzzoha, M. Singh, and D. Lucy. Uncertainty and sensitivity analysis of
the basic reproduction number of a vaccinated epidemic model of influenza. Applied
Mathematical Modelling, 37(3):903–915, 2013.

[21] H. S. Rodrigues, M. T. T. Monteiro, and D. F. M. Torres. Seasonality effects on
dengue: basic reproduction number, sensitivity analysis and optimal control. Math-
ematical Methods in the Applied Sciences, 39(16):4671–4679, 2016.

[22] V. Capasso and G. Serio. A generalization of the kermack-mckendrick deterministic
epidemic model. Mathematical Biosciences, 42(1-2):43–61, 1978.

[23] S. R. Chawla, S. Ahmad, W. Albalawi, A. Khan, I. Shah, and M. R. Eid. Stability
analysis of a modified general seir model with harmonic mean type of incidence rate.
Alexandria Engineering Journal, 127:1183–1192, 2025.

[24] A. Korobeinikov and P. K. Maini. Non-linear incidence and stability of infectious
disease models. Mathematical Medicine and Biology, 22(2):113–128, 2005.