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

A Mathematical Model of COVID-19 Using Piecewise
Derivative of Fractional Order

Shabana Naz1, Muhammad Sarwar1,2, Kamal Shah1,2, Nahid Fatima2,
Thabet Abdeljawad3,2,4,∗

1 Department of Mathematics, University of Malakand, Chakdara Dir(L), 18000,
Khyber Pakhtunkhwa, Pakistan
2 Department of Mathematics and Sciences, Prince Sultan University, Riyadh, 11586,
Saudi Arabia
3 Department of Mathematics, Saveetha School of Engineering, Saveetha Institute of Medical
and Technical Sciences, Saveetha University, Chennai 602105, Tamil Nadu, India
4 Department of Mathematics and Applied Mathematics, School of Science and
Technology, Sefako Makgatho Health Sciences University, Ga-Rankuwa, South Africa

Abstract. Currently the dynamical systems of infectious disease were studied by using various
definitions of fractional calculus. Because the mentioned area has the ability to demonstrate the
short and long memory terms involved in the physical dynamics of numerous real world problems.
In this work, we consider a seven compartmental model for the transmission dynamics of COVID-19
including susceptible (S), vaccinated (V ), exposed (E), infected (I), quarantined (Q), recovered
(R), and death (D) classes. We first revisit the fundamental outcomes related to equilibrium
points, basic reproduction numbers, sensitivity analysis, and equilibrium points including disease
free (DF) and endemic equilibrium (EE) points. In addition, our main goal is to investigate the
considered model under the new aspect of fractional calculus known as piecewise fractional order
operators. The mentioned operators have the ability to demonstrate the multi phase behaviours
of the dynamical problems. The said characteristics cannot be described by using the traditional
operators. We apply the tools of nonlinear analysis to deduce sufficient conditions for the existence
of at least one solution and its uniqueness. Additionally, we also investigate the results related
to stability analysis of Ulam-Hyers (UH) type. Finally, we extend the concepts of RK2 method
to form a sophisticated algorithm to simulate the results graphically. We present the numerical
results for various fractional orders.

2020 Mathematics Subject Classifications: 26A33, 34A08, 03C65

Key Words and Phrases: Crossover effect, Fractional calculus, Stability, Numerical Results,
Sensitivity

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.5861

Email addresses: sarwarswati@gmail.com (M.Sarwar),
kamalshah408@gmail.com kshah@psu.edu.sa (K. Shah),
nfatima@psu.edu.sa (N. Fatima),tabdeljawad@psu.edu.sa (T. Abdeljawad)

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



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 2 of 33

1. Introduction

Mathematical models constitutes an interesting and applicable area of research to de-
scribe various real world problems. The said concept was initiated Bernoulli in 1776.
Later on the idea was formally presented by Karmark and his co-authors in 1927 who
established a simple susceptible infected and recovered (SIR) type model [13]. In sub-
sequent years the mentioned area was increasingly extended to model various real world
problems mathematically [10]. It is authenticating that epidemiology is the important
area of medical field. Therefore, researchers have increasingly used mathematical models
of classical, or difference type differential equations for mathematical models of various
infectious diseases [11, 16]. Since infectious diseases have been investigated by using the
idea of mathematical models. Researchers in majority of published work have used ordi-
nary derivatives or simple algebraic type equations [19] to study the mentioned area. One
of the major infectious diseases is caused by coronavirus known as severe acute respira-
tory syndrome coronavirus 2 ( SARs-COV-2). Recently the said infectious disease gave
birth a world wide pandemic in which more than 6 millions people have died. Also more
than 60 millions have infected . Still the disease is in progress around the globe [44] . In
the end of 2019, the first patient was reported in Wuhan China. Then many cases were
reported in the coming few months. At the end of March there were 11 millions people
been infected in the Wuhan city in which nearly 4 thousands were died. In April 2020 the
said diseases were spread in many countries of the world and hence WHO announced it
a global outbreak [47]. Many countries of the globe implemented lockdown and imposed
strict precautionary measures to control the disease from further spreading. In the mean
time, researchers of medical and bio-engineering fields focused on investigating the cure
or proper medicine for the diseases [6].

Here, we remark that researchers of bio-engineering and applied sciences developed
various mathematical models to investigate the transmission dynamics of the mentioned
diseases [18]. Researchers used classical calculus and algebraic tools to investigate the said
infectious disease. They used traditional tools for numerical simulations in their models.
For instance authors [29] studied the following compartmental model of SARS-COV-2



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 3 of 33

using integer order derivative as follows:

d

dt
[S(t)] = Λ−

[
b(1− pξ)I(t) + (η + d)

]
S(t),

d

dt
[V (t)] = ηS(t)−

[
wbI(t) + d

]
V (t),

d

dt
[E(t)] = wbV (t)I(t) + b(1− pξ)S(t)I(t)− (τ + d)E(t),

d

dt
[I(t)] = τE(t)− (αI + δI + σI + d)I(t),

d

dt
[Q(t)] = αII(t)− [σQ + δQ + d]Q(t),

d

dt
[R(t)] = σII(t) + σQQ(t)− dR(t),

d

dt
[D(t)] = δII(t) + δQQ(t),

S(0) = S0, V (0) = V0, E(0) = E0, I(0) = I0, Q(0) = Q0,

R(0) = R0, D(0) = D0.

(1)

The nomenclatures are described as: Λ is denoted birth or recruitment rate, b stands for
transmission rate of infection, p for rate of individuals who wear face masks, ξ is used for
the rate of effectiveness of face masks, αI stands for the rate of isolation of infected people.
In addition, σI denoted recovered people rate, δI denotes mortality rate of infection, σQ
denotes recovered people during isolation, δQ represents death rate of quarantined people
due to infection. Further, d stands for natural death rate, and w stands for the rate
of reduction due to vaccination process, η denotes rate of vaccination, and τ stands for
incubation period.

On the other hand due to the significant use of fractional calculus in mathematical
models of real world problems, researchers have introduced many operators for integra-
tions and differentiations with real and complex orders. Since traditional tools of calculus
cannot describe the past history and memory terms involved in the mathematical models
of real world problems. Also, most of the biological populations dynamical systems consist
of discrete type quantities which cannot be understood using ordinary derivatives and inte-
grals. Therefore, researchers recommend for better understanding the dynamical behavior
of real world problems by using the techniques of fractional calculus [15]. Some valuable
work recently conducted on the aforementioned area, we refer to [5, 25, 26, 32, 37]. For
the importance of mathematical models, we refer to [28].

After looking at the existing literature, we found that numerous mathematical models
have been developed to examine the illness of SARS-COV-2 using different concepts of
fractional orders derivatives. As a result, a great deal of research has been done on these
pandemic using compartmental models with vaccinated class. However, after a careful
examination, we noted that the mathematical models of the said disease with seven com-
partments have not been adequately examined in relation to the piecewise derivatives to



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 4 of 33

investigate the multi-phase dynamics. Therefore, in this regard, we need further investi-
gation of such models to develop the computational as well as theoretical results for better
understanding. As of October 10, 2020, the novel coronavirus, known as COVID-19, had
spread to more than 200 nations globally and resulted in over 36 million confirmed cases.
As a result, a number of machine learning models that can predict the outbreak world-
wide have been made public. Researchers have developed huge amount of research work
about the said diseases by using the concepts of mathematical models. Compartmental
models have played great roles in predicting the future dynamics of such infectious dis-
eases. Infectious diseases were regularly investigated by researchers by using various tools
of mathematics. We refer some remarkable work as [22, 31, 35, 41, 48]. The mentioned
area has been greatly applied in various disciplines to study various real world problems.
Because it has many important applications. We refer to some [23, 27, 42, 48]. But frac-
tional derivatives have not been uniquely defined yet. There are several approaches which
are used by researchers. Recently, researchers found that some evolution processes suffer
from abrupt changes in their state of dynamics. In such effect the process starts to show
multi phase behaviours which cannot be explained by traditional calculus or fractional
calculus concepts. Therefore, Atangana and his co-author [8] used piecewise concepts
and brilliantly deduced some remarkable results for various dynamical problems. Some
contribution of piecewise derivatives, we offer here as [20, 21, 45].

Many real world problems exhibit abrupt changes in their state of dynamics. For ex-
ample earth quake, fluctuation in economy of less developed countries, and propagations
of infections in society [46]. The mentioned sudden or abrupt changes in the dynamics of
real world problems cannot be accommodate by using the classical tools of calculus. For
these purposes recently researchers [7] introduced the concept of piecewise derivatives of
fractional orders. The mentioned concepts were applied in many real world problems in-
vestigations using mathematical models. Researchers found that the said area excellently
accommodate the concerned gap and comprehensively explain the dynamics of those evo-
lution processes preserving the crossover or multi-phase dynamics [36]. The mentioned
tools have not been yet applied in the mathematical model of SARS-COV-2 given in (1)
to investigate the crossover dynamics. To fill this gap and to investigate the crossover
dynamics of the model under the piecewise derivative concepts, we will follow the theory
given in [40] and numerical procedures to achieve our study. Also, for numerical tools and
qualitative analysis we utilized the procedures given in [9, 30, 38].

Hence inspired from the useful applications of piecewise derivatives of fractional order,
we will extend the seven compartmental model given in (1) under the piecewise derivative
of fractional order for further study in this project. We will establish the qualitative
theory of existence of solution by using fixed point theory of Banach and Schauder type.
Additionally, we will study the global stability under the mentioned piecewise derivative
by using Volterra-Laypunov method. We also establish a numerical scheme based on the
modified Euler method to simulate our theoretical results graphically. Here, we state that
the local stability around the best numerical solution for any solution of the considered
model (1) will be studied by using UH concept. To achieve the above statements, our



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 5 of 33

proposed model will be described as follows:

PC
0 Dε

t [S(t)] = Λ−
[
b(1− pξ)I(t) + (η + d)

]
S(t),

PC
0 Dε

t [V (t)] = ηS(t)−
[
wbI(t) + d

]
V (t),

PC
0 Dε

t [E(t)] = wbV (t)I(t) + b(1− pξ)S(t)I(t)− (τ + d)E(t),
PC
0 Dε

t [I(t)] = τE(t)− (αI + δI + σI + d)I(t),
PC
0 Dε

t [Q(t)] = αII(t)− [σQ + δQ + d]Q(t),
PC
0 Dε

t [R(t)] = σII(t) + σQQ(t)− dR(t),
PC
0 Dε

t [D(t)] = δII(t) + δQQ(t),

(2)

where PC
0 Dε

t denotes piecewise derivative with fractional order ε ∈ (0, 1]. Here, we in-
vestigate the model over [0,T] which can be expressed in sub intervals described by
∇1 = [0, t1], ∇2 = (t1,T]. The flow chart of the model is given in figure 1.

Figure 1: Flow chart of model (2).

In addition, the involve parameters have been described with their numerical values
as in table 1



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 6 of 33

Nomenclatures Description Values

Λ Natural birth rate 0.0020
b Infection transmission rate 0.00005
p The proportional of facing mask 0.01
ξ The rate of efficacy of face masks 0.01
η Rate of vaccination 0.02202643
τ Period of incubation 0.05
d Natural martially rate 0.019
ω Recovery rate 0.0099259
αI Isolation rate of infected people 0.020
δI The rate at which infected die 0.010
σI The rate at which infected recovered 0.10
σQ The rate at which quarantine recovered 0.013978
δQ The rate at which quarantine die 0.015

Table 1: Nomenclature, description and values.

In our study, we also compute the feasible region, positivity and boundedness of solu-
tion under the mentioned derivative and equilibrium points. Also, the basic reproductive
numbers together with the equilibrium points will be deduced for the model (2). Also
an attempt will be done to investigate the sensitivity analysis of the basic reproduction
number by following [34]. Results related to global stability of endemic equilibrium have
been studied by using Volterra-Lyapunov method [43]. Some results devoted to existence
theory and Ulam stability have been investigated by using tools used in [1, 14, 17] and
[12, 33, 39]. The numerical investigations are investigated by establishing a numerical
scheme based on the procedure given in [24]. Finally the results for various fractional
orders will be presented graphically. Various models of dynamical systems like studied in
[2–4], researchers have used fractional derivatives with power law type kernels to investi-
gate the numerical solutions. Our study uses piecewise version of fractional derivative to
investigate the existence and numerical solution for seven compartmental models.

2. Fundamental Results

Here, we recall some basic definitions from fractional calculus. For the given definitions
we have used [32].

Definition 1. The ML function is a recognisee generalization of the exponential. The ML
function with one argument is defined as follows:

Eε(t) =

∞∑
k=0

tk

Γ(kε+ 1)
,



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 7 of 33

while for double parameters it is defined as described by:

Ep,q(t) =
∞∑
k=0

tk

Γ(kp+ q)
.

Definition 2. [5] The Riemann-Liouville (RL) fractional integral with fractional order
ε > 0 for a continuous function ϕ on [0,T] is defined by

Iε0+ϕ(t) =
1

Γ(ε)

∫ t

0
(t− φ)ε−1ϕ(φ)dφ, (3)

such that integral exists on right side.

Definition 3. [5] The RL differentiation with non-integer order n − 1 < ε < n of a
function ϕ which has absolutely continuous derivative up (n− 1) is defined by

Dε
0+ϕ(t) =

1

Γ(n− ε)

(
d

dt

)n ∫ t

0
(t− φ)n−ε−1ϕ(φ)dφ,

where n = [ε] + 1.

Definition 4. [26] The Caputo fractional order derivative of order n − 1 < ε < n of a
differentiable function ϕ ∈ Cn[0,T] is defined by

CDε
0+ϕ(t) =

1

Γ(n− ε)

∫ t

0
(t− φ)n−ε−1ϕ(n)(φ)dφ,

where n = [ε] + 1.

Lemma 1. [26]The given equation

CDε
0+ϕ(t) = χ(t),

ϕ(0) = ϕ0

(4)

has a solution of the form

ϕ(t) = ϕ0 +
1

Γ(ε)

∫ t

0
(t− φ)ε−1χ(φ)dφ. (5)

Definition 5. The piecewise integration of a continuous function ϕ such that 0 < t1 <
t2 < t3 < t is defined by ∫ t3

t1

ϕ(t)dt =

∫ t2

t1

ϕ(t)dt+

∫ t3

t2

ϕ(t)dt.

Definition 6. [37] The piecewise RL fractional integral of order ε ∈ (0, 1] of a continuous
function ϕ is described as

PRLIε0+ϕ(t) =


∫ t1

0
ϕ(t)dt, with t ∈ ∇1,

1

Γ(ε)

∫ t

t1

(t− φ)ε−1ϕ(φ)dφ, with t ∈ ∇2.



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 8 of 33

Definition 7. [37] The piece-wise derivative of order ε ∈ (0, 1] in the Caputo sense of a
differentiable function ϕ ∈ C[0, T ] is defined by

PC
0 Dε

tϕ(t) =


dϕ(t)

dt
, with t ∈ ∇1,

CDε
0+ϕ(t), with t ∈ ∇2,

such that PCDε
0+ represents ordinary differentiation in ∇1 and the Caputo differentiation

on ∇2.

Lemma 2. [25] The solution of the problem

PC
0 Dε

tϕ(t) = χ(t),

ϕ(0) = ϕ0,

is given by

ϕ(t) =


ϕ0 +

∫ t1

0
χ(φ)dφ, with t ∈ ∇1,

ϕ(t1) +
1

Γ(ε)

∫ t

t1

(t− φ)ε−1χ(φ)dφ, with t ∈ ∇2.

Theorem 1. [14]Let D ⊂ B be closed convex bounded set such that D = {z ∈ B z =
λNz, λ ∈ [0, 1]},, then N : D → B is completely continuous operator has at least one
fixed point.

Definition 8. According to Ulam [12, 33, 39] ”near every exact function there are many
approximation of it”. According to their definition consider a functional equation given by

Nz = z

is UH stable if for any ε > 0, for the inequality

|z −Nz| ≤ ε, z ∈ B,

and constant LN > 0 one has
|z − z∗| ≤ LNε,

where z∗ is the best approximate or exact solution.

Definition 9. Sensitivity index[34]
It is an important threshold value calculated by the formula given by

SR0
p =

p

R0

∂R0

∂p
, (6)

were p is the parameter describing the relation of R0. Sensitivity index tells us how much
be the dynamic effected by increasing or decreasing the values of parameters.



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 9 of 33

3. Fundamental Results of Model (2)

Here we provide the positivity, boundedness, DF and EE points. Also we compute the
basic reproductive number and presents its 3D profile. Sensitivity is also discussed here.

Theorem 2. All solutions of model (2) are bounded and inclined to a feasible region
described by

Ω =

{
(S, V,E, I,Q,R) ∈ R6

+S + V + E + I +Q+R ≤ Λ

d

}
.

Proof. Let M be the total pollution at time t > 0, such that

M(t) = S(t) + V (t) + E(t) + I(t) +Q(t) +R(t). (7)

Applying PC
0 Dε

t to both sides of (7) gives

PC
0 Dε

t [M(t)] = PC
0 Dε

t [S(t)] +
PC
0 Dε

t [V (t)] +PC
0 Dε

t [E(t)] +PC
0 Dε

t [I(t)] (8)

+ PC
0 Dε

t [Q(t)] +PC
0 Dε

t [R(t)].

Using the corresponding values from (2) in (7) and simplifying, we get

PC
0 Dε

t [M(t)] = Λ− d(S + E + V + I +Q+R)

= Λ− dM. (9)

We can write (9)

PC
0 Dε

t [M(t)] =
dM
dt = Λ− dM, t ∈ ∇1,

C
0D

ε
t [M(t)] = Λ− dM, t ∈ ∇2.

(10)

Applying Laplace transform in (15) yields that

M(t) =
Λ
d −M0 exp(−dt), t ∈ ∇1,

Λ
d (1−Eε(−dtε)) +Eε(−dtε), t ∈ ∇2.

(11)

Hence at t → ∞, we get from (11)

M(t) ≤ Λ

d
.

Theorem 3. All the solutions (S, V,E, I,Q,R,D) of model (2) are non-negative for t > 0
with positive initial condition (S0, V0, E0, I0, Q0, R0, D0) > 0.



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 10 of 33

Proof. From model (2), we have

PC
0 Dε

t [S(t)]

∣∣∣∣
S=0

= Λ > 0,

PC
0 Dε

t [V (t)]

∣∣∣∣
V=0

= 0,

PC
0 Dε

t [E(t)]

∣∣∣∣
E=0

= 0,

PC
0 Dε

t [I(t)]

∣∣∣∣
I=0

= 0,

PC
0 Dε

t [Q(t)]

∣∣∣∣
Q=0

= 0,

PC
0 Dε

t [R(t)]

∣∣∣∣
R=0

= 0,

PC
0 Dε

t [D(t)]

∣∣∣∣
D=0

= 0,

from (12), one see that PC
0 Dε

t [S(t)] > 0 which yields that S(t) > 0, at t > 0, hence we
conclude that (S, V,E, I,Q,R,D) > 0 at t > 0. Hence positivity condition derived.

Further, the DFE and endemic equilibrium points can be computed by using the results
of [29] as follows:
DFE
From (2) putting left sides equal to zeros, that is

Λ−
[
b(1− pξ)I(t) + (η + d)

]
S(t) = 0,

ηS(t)−
[
wbI(t) + d

]
V (t) = 0,

wbV (t)I(t) + b(1− pξ)S(t)I(t)− (τ + d)E(t) = 0,

τE(t)− (αI + δI + σI + d)I(t) = 0,

αII(t)− [σQ + δQ + d]Q(t) = 0,

σII(t) + σQQ(t)− dR(t) = 0,

δII(t) + δQQ(t) = 0.


Now I = 0, E = 0, Q = 0, R = 0, D = 0 due to absence of infection, and S = S0, V = V 0,
and solving system (12), we have

E0 =

(
Λ

η + d
,

Λ

d(η + d)
, 0, 0, 0, 0, 0

)
.

Endemic equilibrium point

In the same way, we can compute endemic equilibrium point E∗ =
(
S̃, Ṽ , Ẽ, Ĩ, Q̃, R̃, D̃

)
,



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 11 of 33

again from equating right sides of model (2) equal to zero as follows:

Λ−
[
b(1− pξ)Ĩ(t) + (η + d)

]
S̃(t) = 0,

ηS̃(t)−
[
wbĨ(t) + d

]
Ṽ (t) = 0,

wbṼ (t)Ĩ(t) + b(1− pξ)S̃(t)Ĩ(t)− (τ + d)Ẽ(t) = 0,

τ Ẽ(t)− (αI + δI + σI + d)Ĩ(t) = 0,

αI Ĩ(t)− [σQ + δQ + d]Q̃(t) = 0,

σI Ĩ(t) + σQQ̃(t)− dR̃(t) = 0,

δI Ĩ(t) + δQQ̃(t) = 0.


On solving (12), we get

S̃ =
Λ

b(1− pξ)Ĩ + η + d
,

Ṽ =
η

wbĨ + d

[
Λ

b(1− pξ)Ĩ + η + d

]
,

Ẽ =
1

θ + d

[
Λb(1− pξ)

b(1− pξ)Ĩ + η + d
+

wbηΛ

b(1− pξ)Ĩ + η + d

]
Ĩ ,

Ĩ =
θ

(αI + δI + σI + d)

[
1

θ + d

[
Λb(1− pξ)

b(1− pξ)Ĩ + η + d
+

wbηΛ

b(1− pξ)Ĩ + η + d

]
Ĩ

]
,

Q̃ =
αI Ĩ

σQ + δQ + d
,

R̃ =
σI
d
Ĩ +

σQ
d

αI Ĩ

σQ + δQ + d
,

Q̃ = 0.

The basic reproductive number R0 has been computed for integer order model in [29], we
recall it as follows:

Rvac =
τb(1− pξ)Λ

(τ + d)(αI + δI + σI + d)

[
d+ ω

d(η + d)

]
, (12)

but if no vaccination, then ω = 0, we get the fundamental reproductive number as follows:

R0 =
τb(1− pξ)Λ

d(τ + d)(αI + δI + σI + d)
. (13)

In figure 2, we describe surface plot of the reproductive numbers to see the behaviour on
decreasing or increasing the parameters values.



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 12 of 33

10.80.6

Λ

0.40.200
0.2

0.4
ξ

0.6
0.8

0

0.02

0.01

0.03

1

R
0

1
0.8

0.6

d
0.4

0.2
00

0.5

ξ

0

1

0.5

1

×10-5

R
0

1
0.8

0.6

Λ

0.4
0.2

00
0.2

0.4
τ

0.6
0.8

0

0.01

0.02

0.03

0.04

1

R
0

10.80.6

p
0.40.200

0.5
τ

8

7

5

6

1

×10-5

R
0

Figure 2: Surface plots of R0.

We see that as the values of the given parameters are decreasing, the value of R0 is
also increasing but never exceeds 1.

3.1. Sensitivity Analysis

Next, we apply the formula to get the sensitivity of the fundamental reproduction
number.

SR0
q =

q

R0

∂R0

∂q
, (14)

where q represents parameter involve in the R0. According to (14), we compute sensitivity
indices 2.

Sensitivity index Value Sensitivity index Value

SR0
b 1.0 SR0

τ 1.0

SR0
ξ - 0.637 SR0

d -0.4567

SR0
p -0.776 SR0

αI
-0.7675

SR0
Λ 1.0 SR0

δI
-0.9986

SR0
σI

-0.895

Table 2: Sensitivity indices of basic reproduction number.

From sensitivity indices table 2, we see that increase of specific percentage of values
means that percentage increase or decrease in the index which implies the impact of
parameters values for simulations in the current studies. Here, in figure 3, we present
the pi chart of sensitivity indices which describe the percentage impact of variable on the
reproductive number.



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 13 of 33

Figure 3: Plot of Sensitivity index.

4. Existence Theory

Let J = [0,T], 0 ≤ t ≤ T < ∞, and ϕ = (S, V,E, I,Q,R,D), then the right sides of
(2), we can express as

L1(t, ϕ(t)) = Λ−
[
b(1− pξ)I(t) + (η + d)

]
S(t),

L2(t, ϕ(t)) = ηS(t)−
[
wbI(t) + d

]
V (t),

L3(t, ϕ(t)) = wbV (t)I(t) + b(1− pξ)S(t)I(t)− (τ + d)E(t),

L4(t, ϕ(t)) = τE(t)− (εI + δI + σI + d)I(t),

L5(t, ϕ(t)) = εII(t)− [σQ + δQ + d]Q(t),

L6(t, ϕ(t)) = σII(t) + σQQ(t)− dR(t),

L7(t, ϕ(t)) = δII(t) + δQQ(t).

(15)

Utilizing (15), model (2), with ε ∈ (0, 1] may be express as:{
PC
0 Dε

tϕ(t) = g(t, ϕ(t)), t ∈ J ,

ϕ(0) = ϕ0,
(16)

where

ϕ(t) =



S(t),

V (t),

E(t),

I(t),

Q(t),

R(t),

D(t).

, ϕ0(t) =



S0,

V0,

E0,

I0,

Q0,

R0,

D0.

g(t, ϕ(t)) =



L1(t, ϕ(t)),

L2(t, ϕ(t)),

L3(t, ϕ(t)),

L4(t, ϕ(t)),

L5(t, ϕ(t)),

L6(t, ϕ(t)).

(17)



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 14 of 33

Lemma 3. For the function g : [0,T×R → R, the given problem

PC
0 Dε

tϕ(t) = g(t, ϕ(t)), 0 < ε ≤ 1, (18)

ϕ(0) = ϕ0

is equivalent to the integral equation described by

ϕ(t) =


ϕ0 +

∫ t

0
g(u, ϕ(u))du, t ∈ ∇1,

ϕ(t1) +

∫ T

t1

(t− u)ε−1

Γ(ε)
g(u, ϕ(u))du, t1 < t ≤ T,

(19)

Proof. Following the fashion of Lemma 2, we can easily obtain the solution above.

The hypothesis given hold:

(H1) Let Lg > 0 is real number and ϕ, ϕ̄ ∈ B, then

|g(t, ϕ)− g(t, ϕ̄)| ≤ Lg|ϕ− ϕ̄|.

(H2) Let Cg > 0 and Mg > 0, then

|g(t, ϕ(t))| ≤ Cg|ϕ|+Mg.

Theorem 4. At least one solution exists for the Problem (18) under the hypotheses
(H1), (H2).

Proof. Let define a closed and convex subset D ⊂ B as

D = {ϕ ∈ B : ∥ϕ∥ ≤ r1,2, r1,2 > 0}.

Let N : D → D using (35), then

N(ϕ) =


ϕ0 +

∫ t

0
g(u, ϕ(u))du, t ∈ ∇1,

ϕ(t1) +

∫ T

t1

(t− u)ε−1

Γ(ε)
g(u, ϕ(u))du, t1 < t ≤ T.

(20)

For ϕ ∈ D, it follows that

|N(ϕ)(t)| ≤


|ϕ0|+

∫ t1

0
|g(u, ϕ(u))|du,

|ϕ(t1)|+
∫ T

t1

(t− u)ε−1

Γ(ε)
|g(u, ϕ(u))| du



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 15 of 33

≤


|ϕ0|+

∫ t1

0
[Cg∥ϕ∥+Mg] du,

|ϕ(t1)|+
∫ T

t1

(t− u)ε−1

Γ(ε)
[Cg∥ϕ∥+Mg] du,

≤


|ϕ0|+ t1 [Cgr1,2 +Mg] ≤ r1,2, t ∈ ∇1,

|ϕ(t1)|+
Tε

Γ(ε+ 1)
[Cgr1,2 +Mg] ≤ r1,2, t1 < t ≤ T,

for t1 < t ≤ T, use |(t1 − u)ε − (T− u)ε| ≤ Tε, one has

r1,2 ≥ max


|ϕ0|+ t1Mg

1− t1Cg
, t ∈ ∇1,

|ϕ(t1)|Γ(ε+ 1) + TεMg

(Γ(ε+ 1)− TεCg
, t1 < t ≤ T.

Thus proved that ∥N(ϕ)∥ ≤ r1,2, implies that N(D) ⊂ D. Thus, bounded set to bounded
set is mapped by N. Hence, N is bounded. Take tm < tn ∈ [0, t1], we have

|N(ϕ)(tn)−N(ϕ)(tm)| =

∣∣∣∣ ∫ tn

0
g(u, ϕ(u))du−

∫ tm

0
g(u, ϕ(u))du

∣∣∣∣
≤

∫ tn

tm

|g(u, ϕ(u))|du

≤
∫ tn

tm

[Cg∥ϕ∥+Mg] du

≤ (Cgr1,2 +Mg)[tn − tm]. (21)

The (21) tells us that tm → tn, implies that

|N(ϕ)(tn)−N(ϕ)(tm)| → 0.

The uniform continuity of N implies that

∥N(ϕ)(tn)−N(ϕ)(tm)∥ → 0, as tn → tm.

It follows that N is equi-continuous in this case.
Next, we take the other interval tm < tn ∈ (t1,T] as

|N(ϕ)(tn)−N(ϕ)(tm)| =

∣∣∣∣ ∫ tn

0

(tn − u)ε−1

Γ(ε)
g(u, g(u))du−

∫ tm

0

(tm − τ)ε−1

Γ(ε)
g(u, ϕ(u))du

∣∣∣∣
≤

∫ tm

0

[(tm − u)ε−1 − (tn − u)ε−1]

Γ(ε)
|g(u, ϕ(u))|du

+

∫ tn

tm

(tn − u)ε−1

Γ(ε)
|g(u, ϕ(u))|du



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 16 of 33

≤
[ ∫ tm

0

[(tm − u)ε−1 − (tn − u)ε−1]

Γ(ε)
du

+

∫ tn

tm

(tn − u)ε−1

Γ(ε)
du

]
(Cg∥ϕ∥+Mg)

≤ (Cgr1,2 +Mg)

Γ(ε+ 1)
[tεn − tεm + 2(tn − tm)

ε] . (22)

We see from (22), that

|N(ϕ)(tn)−N(ϕ)(tm)| → 0, as tm → tn.

As N is bounded over [t1,T] and also uniformly continuous. Thus

∥N(ϕ)(tn)−N(ϕ)(tm)∥ → 0, as tm → tn.

Hence, N is equi-continuous over (t1,T].
Therefore, N is equi- continuous mapping over [0, t1]∪(t1,T]. By Arzelá-Ascoli theo-

rem, N is completely continuous. Thus in view of Theorem 1, the problem (18) has at
least one solution.

Theorem 5. Under the hypothesis (H1), the problem (18) has a unique solution with

max

{
t1Lg,

Tε

Γ(ε+ 1)
Lg

}
= ∆ < 1

holds.

Proof. Let N : B → B be the mapping, then for ϕ, ϕ̄ ∈ B, we consider for [0, t1] as

∥N(ϕ)−N(ϕ̄)∥ = max
t∈[0,t1]

∣∣∣∣ ∫ t1

0
g(u, ϕ(u))du−

∫ t1

0
g(u, ϕ̄(u))du

∣∣∣∣
≤ t1Lg∥ϕ− ϕ̄∥. (23)

Hence , we have from (23)

∥N(ϕ)−N(ϕ̄)∥ ≤ t1Lg∥ϕ− ϕ̄∥. (24)

If t1 ≤ t ≤ T, one has

∥N(ϕ)−N(ϕ̄)∥ = max
t∈(t1,T]

∣∣∣∣ ∫ T

t1

(t− u)ε−1

Γ(ε)
g(u, ϕ(u))du−

∫ T

t1

(t− u)ε−1

Γ(ε)
g(u, ϕ̄(u))du

∣∣∣∣
≤ TεLg

Γ(ε+ 1)
∥ϕ− ϕ̄∥. (25)

Let max

{
TLg,

TεLg

Γ(ε+1)

}
= ∆. Then, from (24) and (25)), we have

∥N(ϕ)−N(ϕ̄)∥ ≤ ∆∥ϕ− ϕ̄∥. (26)

As a result, an operator of contraction is N. As a result, there is only one solution to the
Problem (18).



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 17 of 33

5. Stability Theory

Here we discuss the stability results for local and global stability as well as UH stability.

5.1. Stability Results

We use the Jacobian matrix method and Lyapunov- Volterra method to prove some
results for LAS and GAS of DFE and endemic equilibrium points.

Theorem 6. The DFE point is LAS if R0 < 1, otherwise unstable.

Proof. The Jacobian matrix say J at DFE
(

Λ
η+d ,

Λ
d(η+d) , 0, 0, 0, 0, 0

)
of model (2) is

given by

J =



−(η + d) 0 0 −b(1− pξ) Λ
η+d 0 0

η −d 0 −wbΛ
d(d+η) 0 0

0 0 −(τ + d) wbΛ
d(d+η) +

b(1−pξ)Λ
η+d 0 0

0 0 τ −(αI + σI + d) 0 0
0 0 0 σI −(σQ + δQ + d) 0
0 0 0 σI σQ −d


(27)

Computing the eigenvalues of (27), we have

(λ+ (η + d))(λ+ µ)(λ4 + a1λ
3 + a2λ

2 + a3λ+ a4) = 0, (28)

where

a1 = k + l + τ + d, a2 = kl +m+ n+ (k + l)(τ + d), a3 = km+ (τ + d)(kl +m) + nl, a4 = km(τ + d) + nm,

such that

k = σI + αI + d, l = σQ + δQ + 2d, m = d(σQ + δQ + d),

n =
(τb(w + a(1− pξ))Λ

d(η + d)
.

We see from (28) that λ1 = −(η + d) < 0, λ2 = −µ < 0, which implies that real parts of
eigenvalues are negative. It is possible to demonstrate that the eigenvalues have negative
real parts by applying the Routh-Hurwitz criterion. Hence, we get from

λ4 + a1λ
3 + a2λ

2 + a3λ+ a4 = 0,

from which we see that

λ4 : a0 a2 a4

λ3 : b1 b3 b5

λ2 : c1 c3 c5



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 18 of 33

λ1 : d1 d3 d5

λ0 : e1 e3 e5,

where

b1 =
−1

a1

∣∣∣∣ a0 a2
a1 a3

∣∣∣∣ = −(a3 − a1a2)

a1
, b3 =

−1

a1

∣∣∣∣ a0 a4
a1 a5

∣∣∣∣ = −(a0a5 − a1a4)

a1
,

c1 =
−1

b1

∣∣∣∣ a1 a3
b1 b3

∣∣∣∣ = −(a1b3 − a3b1)

b1
, c3 =

−1

b1

∣∣∣∣ a1 a5
b1 b5

∣∣∣∣ = −(a1b5 − a5b1)

b1
,

d1 =
−1

c1

∣∣∣∣ a1 a3
c1 c3

∣∣∣∣ = −(a1c3 − a3c1)

c1
, d3 =

−1

c1

∣∣∣∣ a1 a5
c1 c5

∣∣∣∣ = −(a1c5 − a5c1)

c1
,

e1 =
−1

d1

∣∣∣∣ a1 a3
d1 d3

∣∣∣∣ = −(a1d3 − a3d1)

d1
, e3 =

−1

d1

∣∣∣∣ a1 a5
d1 d5

∣∣∣∣ = −(a1d5 − a5d1)

d1
.

By computing the above determinants, we can see that the coefficients in the first column
are of same positive sign. Hence, all the roots have negative real parts. Thus the DFE is
locally asymptotically stable.

Theorem 7. The EE point E∗ of the proposed model (2) is GAS in the given feasible
region Ω.

Proof. Let us define a Lyapunov function F : Ω → R by

F (t) = q1

(
S − S∗ − S∗ ln

[
S

S∗

])
+ q2

(
V − V ∗ − V ∗ ln

[
V

V ∗

])
+ q3

(
E − E∗ − E∗ ln

[
E

E∗

])
= q4

(
I − I∗ − I∗ ln

[
I

I∗

])
+ q5

(
Q−Q∗ −Q∗ ln

[
Q

Q∗

])
(29)

+ q6

(
R−R∗ −R∗ ln

[
R

R∗

])
+ q7

(
D −D∗ −D∗ ln

[
D

D∗

])
,

where qi, i = 1, 2, ..., 7 are positive constants and will be specified latter. Upon application
of PC

0 Dε
t on both sides of (29), and replacing S, V,E, I,Q,R,D by S − S∗, V − V ∗, E −

E∗, I − I∗, Q−Q∗, R−R∗, D −D8, respectively, we have

PC
0 Dε

t [F (t)] = q1

[
S − S∗

S

]
PC
0 Dε

tS(t) + q2

[
V − V ∗

V

]
PC
0 Dε

tV (t)



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 19 of 33

+ q3

[
E − E∗

E

]
PC
0 Dε

tE(t) + q4

[
I − I∗

I

]
PC
0 Dε

tI(t)

+ q5

[
Q−Q∗

Q

]
PC
0 Dε

tQ(t) + q6

[
R−R∗

R

]
PC
0 Dε

tR(t)

+ q7

[
D −D∗

D

]
PC
0 Dε

tD(t). (30)

By choosing qi = 1, for i = 1, 2, ..., 7, and plugging values from (2) in (30), and separating
positive and negative parts, such that

Θ1 = Λ
S − S∗

S
+ η

(V − V ∗)(S − S∗)

V
+ wb

(E − E∗)(V − V ∗)(I − I∗)

E

+ b(1− pξ)
(S − S∗)(I − I∗)(E − E∗)

E
+ τ

(E − E∗)(I − I∗)

I
+ αI

(Q−Q∗)(I − I∗)

Q
(31)

+ σI(R−R∗)
(I − I∗)

R
+ σQ

(Q−Q∗)(R−R∗)

R
+ δI

(D −D∗)(I − I∗)

D

+ δQ
(Q−Q∗)(D −D∗)

D
,

and

Θ2 = b(1− pξ)
(S − S∗)2(I − I∗)

S
+ (η + d)

(S − S∗)2

S
+ wb

(E − E∗)(V − V ∗)2(I − I∗)

V

+
d(V − v∗)2

V
+ (τ + d)

(E − E∗)2

E
+ (αI + δI + σI + d)

(I − I∗)2

I
(32)

+ (σQ + δQ + d)
(Q−Q∗)2

Q
+ d

(R−R∗)2

R
.

Using (31), and (32) in (30), we have

PC
0 Dε

t [F (t)] = Θ1 −Θ2. (33)

We see that in (33) if Θ1 < Θ2, then
PC
0 Dε

t [F (t)] ≤ 0, Hence the Lyapunov- Volterra
function F is negative under the derivative of piecewise. But F > 0, therefore F is
positive definite function, therefore, the endemic equilibrium point E∗ is GAS.

5.2. UH Stability

The UH stability results are derived in this portion.

Remark 1. If h is a function that is independent of ϕ, then h : J → (0,∞), ∋ h(0) =
h(t1) = 0 and

|h(t)| ≤ ϵ, for every t ∈ J .



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 20 of 33

Lemma 4. The

PC
0 Dε

tϕ(t) = g(t, ϕ(t)) + h(t), 0 < ε ≤ 1,

ϕ(0) = ϕ0

yield a result of the form

ϕ(t) =


ϕ0 +

∫ t

0
g(u, ϕ(u))du+

∫ t

0
h(u)du, t ∈ ∇1,

ϕ(t1) +

∫ T

t1

(t− u)ε−1

Γ(ε)
g(u, ϕ(u))d(u) +

∫ T

t1

(t− u)ε−1

Γ(ε)
h(u)du, t ∈ ∇2.

(34)

Moreover, the solution (34) satisfies the following relation in view of Remark 1 and using
|T− t1| ≤ T as

∣∣∣∣ϕ(t)− ϕ0 −
∫ t

0
g(u, ϕ(u))du

∣∣∣∣ ≤ t1ϵ, t ∈ ∇1,∣∣∣∣ϕ(t)− ϕ(t1)−
∫ T

t1

(t− u)ε−1

Γ(ε)
g(u, ϕ(u))du

∣∣∣∣ ≤ Tεϵ

Γ(ε+ 1)
, t ∈ ∇2.

Proof. The proof of the Lemma 4 is easy to derive.

Theorem 8. We get UH and generalized UH stable if the condition max

{
1 − t1Lg, 1 −

∆̄

}
< 1, where ∆̄ =

TεLg

Γ(ε+1) holds corresponding to the solution of problem (18).

Proof. We derive proof in two steps:
Step I: t ∈ [0, t1], then

∥∥ϕ− ϕ̄
∥∥ = max

t∈J

∣∣∣∣ϕ(t)− ϕ0 −
∫ t

0
g(u, ϕ̄(u))du

∣∣∣∣
≤ max

t∈J

∣∣∣∣ϕ(t)− ϕ0 −
∫ t

0
g(u, ϕ(u))du

∣∣∣∣
+ max

t∈J

∣∣∣∣ ∫ t

0
g(u, ϕ(u))du−

∫ t

0
g(u, ϕ̄(u))du

∣∣∣∣
≤ ϵt1 + Lgt1∥ϕ− ϕ̄∥. (35)

Hence from (35), we get ∥∥ϕ− ϕ̄
∥∥ ≤ ϵt1

1− t1Lg
. (36)



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 21 of 33

Step II: If t ∈ [t1,T], then we have

∥∥ϕ− ϕ̄
∥∥ = max

t∈J

∣∣∣∣ϕ(t)− ϕ(t1)−
∫ T

t1

(t− u)ε−1

Γ(ε)
g(u, ϕ̄(u))du

∣∣∣∣
≤ max

t∈J

∣∣∣∣ϕ(t)− ϕ(t1)−
∫ T

t1

(t− u)ε−1

Γ(ε)
g(u, ϕ(u))du

∣∣∣∣
+ max

t∈J

∣∣∣∣ ∫ T

t1

(t− u)ε−1

Γ(ε)
g(u, ϕ(u))du−

∫ T

t1

(t− u)ε−1

Γ(ε)
g(u, ϕ̄(u))du

∣∣∣∣
≤ Tεϵ

Γ(ε+ 1)
+

TεLg

Γ(ε+ 1)
∥ϕ− ϕ̄∥. (37)

From (37), ∥∥ϕ− ϕ̄
∥∥ ≤

[
Tε

(1− ∆̄)Γ(ε+ 1)

]
ϵ. (38)

We can now see that the solution to issue (18) is UH stable based on (36) and (38).
Moreover, the solution is unquestionably generalized UH stable if there is a non-decreasing
function Φ : J → [0,∞) such that Φ(ϵ) = ϵ.

6. Numerical Scheme

Here, we expand the approach from [24] to build the numerical scheme based on
the RKm method of order four. For computation purposes, where small steps are to
be computed the Euler method is good. However, for evaluations of largen numbers of
steps, the RKm method is given preferred because in such a situation it gives efficient and
accurate solution. In light of this need, we view the fractional order problem as

PC
0 Dε

tϕ(t) = gj(t, ϕ(t)), j = 1, 2, ..., 7

ϕ(0) = ϕ0. (39)

where gj : J × R → R. We may express Taylor series in general form as

ϕ(t+ h) = ϕ(t) +
hε

Γ(ε+ 1)
Dt

εϕ(t) +
h2ε

Γ(2ε+ 1)
D2ε

t x(t) + ... (40)

On applying Dt
2εϕ = Dtgj(t, ϕ(t)) + gj(t, ϕ(t))Dϕ

εg(t, ϕ(t)) in (40) implies the following
results:

ϕi+1 = ϕi +
hε

Γ(ε+ 1)
gj(ti, ϕ(ti)) +

hε

2Γ(ε+ 1)

[
K1 +K2

]
, (41)

where

K1 = gj(ti, ϕi(ti)), K2 = gj

(
ti +

2hεΓ(ε+ 1)

Γ(2ε+ 1)
, ϕ(ti) +

2hεΓ(ε+ 1)

Γ(2ε+ 1)
gj(ti, ϕi(ti))

)
. (42)



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 22 of 33

On using ε = 1, gives ordinary numerical method.
To illustrate the compartmental solutions graphically, we derive numerical scheme as:

S(ti+1)) =


Si−1(ti−1) +

h

2
g1

[
ti−1 +

h

2
, ϕi−1(ti−1) +

K1

2

]
, t ∈ ∇1

Si(ti) +
hφ

Γ(ε+ 1)
g1(ti, ϕi(ti)) +

hε

2Γ(ε+ 1)

[
K2 +K3

]
, t ∈ ∇2,

(43)

where h = ti+1 − ti and

K1 = g1(ti−1, ϕi−1(ti−1)), K2 = g1(ti, ϕi(ti)),

K3 = g1

(
ti +

2hεΓ(ε+ 1)

Γ(2ε+ 1)
, ϕ(ti) +

2hεΓ(ε+ 1)

Γ(2ε+ 1)
g1(ti, ϕi(ti))

)
. (44)

With the same patron

V(ti+1)) =


Vi−1(ti−1) +

h

2
g2

[
ti−1 +

h

2
, ϕi−1(ti−1) +

K1

2

]
, t ∈ ∇1

Vi(ti) +
hε

Γ(ε+ 1)
g2(ti, ϕi(ti)) +

hε

2Γ(ε+ 1)

[
K2 +K3

]
, t ∈ ∇2,

(45)

E(ti+1)) =


Ei−1(ti−1) +

h

2
g3

[
ti−1 +

h

2
, ϕi−1(ti−1) +

K1

2

]
, t ∈ ∇1

Ei(ti) +
hε

Γ(ε+ 1)
g3(ti, ϕi(ti)) +

hε

2Γ(ε+ 1)

[
K2 +K3

]
, t ∈ ∇2,

(46)

I(ti+1)) =


Ii−1(ti−1) +

h

2
g4

[
ti−1 +

h

2
, ϕi−1(ti−1) +

K1

2

]
, t ∈ ∇1

Ii(ti) +
hε

Γ(ε+ 1)
g4(ti, ϕi(ti)) +

hε

2Γ(ε+ 1)

[
K2 +K3

]
, t ∈ ∇2,

(47)

Q(ti+1)) =


Qi−1(ti−1) +

h

2
g5

[
ti−1 +

h

2
, ϕi−1(ti−1) +

K1

2

]
, t ∈ ∇1

Ri(ti) +
hε

Γ(ε+ 1)
g5(ti, ϕi(ti)) +

hε

2Γ(ε+ 1)

[
K2 +K3

]
, t ∈ ∇2,

(48)

R(ti+1)) =


Ri−1(ti−1) +

h

2
g6

[
ti−1 +

h

2
, ϕi−1(ti−1) +

K1

2

]
, t ∈ ∇1

Ri(ti) +
hε

Γ(ε+ 1)
g6(ti, ϕi(ti)) +

hε

2Γ(ε+ 1)

[
K2 +K3

]
, t ∈ ∇2,

(49)

and

D(ti+1)) =


Di−1(ti−1) +

h

2
g7

[
ti−1 +

h

2
, ϕi−1(ti−1) +

K1

2

]
, t ∈ ∇1

Di(ti) +
hε

Γ(ε+ 1)
g7(ti, ϕi(ti)) +

hε

2Γ(ε+ 1)

[
K2 +K3

]
, t ∈ ∇2.

(50)



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 23 of 33

7. Graphical Illustrations and Explanations

To simulate our results graphically, we use the numerical values of table 1 and taking
the initial data as given by S(0) = 97.00000; V (0) = 50.00000; E(0) = 70.00000; I(0) =
1.581936; Q(0) = 1.581936; R(0) = 1.538689; D(0) = 0.030664.

7.1. Case-I

We take the first set of fractional order values in (0.60, 0.80] and presented solutions
graphically in figures 4, 5, 6, 7, 8, 9, 10 respectively.

t
0 50 100 150 200 250 300

S
us

ce
pt

ib
le

 c
la

ss

0

20

40

60

80

100
0.65
0.70
0.75
0.80

Figure 4: Solution illustration of susceptible class for ε ∈ (0.60, 0.80].

t
0 50 100 150 200 250 300

V
ac

ci
na

te
d 

cl
as

s

50

60

70

80

90

0.65
0.70
0.75
0.80

Figure 5: Solution illustration of vaccinated class for ε ∈ (0.60, 0.80].



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 24 of 33

t
0 50 100 150 200 250 300

E
xp

os
ed

 c
la

ss
 

0

20

40

60

80
0.65
0.70
0.75
0.80

Figure 6: Solution illustration of exposed class for ε ∈ (0.60, 0.80].

t
0 50 100 150 200 250 300

In
fe

ct
ed

 c
la

ss

0

5

10

15

20

25
0.65
0.70
0.75
0.80

Figure 7: Solution illustration of infected class for ε ∈ (0.60, 0.80].



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 25 of 33

0 50 100 150 200 250 300

t

0

1

2

3

4

5

6

Q
ua

ra
nt

in
ed

 c
la

ss

0.65
0.70
0.75
0.80

Figure 8: Solution illustration of quarantined class for ε ∈ (0.60, 0.80].

t
0 50 100 150 200 250 300

R
ec

ov
er

ed
 c

la
ss

0

2

4

6

8
0.65
0.70
0.75
0.80

Figure 9: Solution illustration of recovered class for ε ∈ (0.60, 0.80].



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 26 of 33

t
0 50 100 150 200 250 300

D
ea

th
 c

la
ss

0

5

10

15

20

25

0.65
0.70
0.75
0.80

Figure 10: Solution illustration of death class for ε ∈ (0.60, 0.80].

In the above graphical presentations, we have taken t1 < 50, and T = 300, the cor-
responding crossover effect can be observed near the point t1 < 50 in all compartments.
The susceptible population is decreasing, in the same vaccinated classes has raised. On
the other hand exposed class was reducing as infection was spreading in the society. In
the mean time quarantined initially risen and then went on on decreasing. The recovered
individuals were increasing with the passage of time and also the death class has grown.
In all compartment stability after some time occurred.

7.2. Case-II

We take the first set of fractional order values in (0.80, 1.00] and presented solutions
graphically in figures 11, 12, 13, 14, 15, 16, 17 respectively.

t
0 50 100 150 200 250 300

S
us

ce
pt

ib
le

 c
la

ss

0

20

40

60

80

100
0.85
0.90
0.95
1.00

Figure 11: Solution illustration of susceptible class for ε ∈ (0.80, 1.00].



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 27 of 33

t
0 50 100 150 200 250 300

V
ac

ci
na

te
d 

cl
as

s

50

60

70

80

90

0.85
0.90
0.95
1.00

Figure 12: Solution illustration of vaccinated class for ε ∈ (0.80, 1.00].

t
0 50 100 150 200 250 300

E
xp

os
ed

 c
la

ss
 

0

20

40

60

80
0.85
0.90
0.95
1.00

Figure 13: Solution illustration of exposed class for ε ∈ (0.80, 1.00].



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 28 of 33

t
0 50 100 150 200 250 300

In
fe

ct
ed

 c
la

ss

0

5

10

15

20

25
0.85
0.90
0.95
1.00

Figure 14: Solution illustration of infected class for ε ∈ (0.80, 1.00].

0 50 100 150 200 250 300

t

0

1

2

3

4

5

6

Q
ua

ra
nt

in
ed

 c
la

ss

0.85
0.90
0.95
1.00

Figure 15: Solution illustration of quarantined class for ε ∈ (0.80, 1.00].



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 29 of 33

t
0 50 100 150 200 250 300

R
ec

ov
er

ed
 c

la
ss

0

2

4

6

8
0.85
0.90
0.95
1.00

Figure 16: Solution illustration of recovered class for ε ∈ (0.80, 1.00].

t
0 50 100 150 200 250 300

D
ea

th
 c

la
ss

0

5

10

15

20

25

0.85
0.90
0.95
1.00

Figure 17: Solution illustration of death class for ε ∈ (0.80, 1.00].

Given t1 < 50 and T = 300 in the graphical displays above, the relevant crossover effect
is visible in all compartments close to t1 < 50. In the same classes that have increased
in vaccinations, the vulnerable population is declining. On the other hand, as infection
increased throughout society, the exposed class shrank. Quarantined cases increased in the
interim before continuing to decline. Over time, both the number of recovered individuals
and the number of deaths have increased. After a while, stability happened in every
compartment. From the figures given above, we see that the decay dynamics is faster at
lower values of fractional orders while in case of growth, the dynamic is rapidly increasing
under the larger values of fractional orders. In the same way, the growth in dynamics is
slower at small fractional orders and faster at higher values. This chastities demonstrate
the influence of fractional orders on the dynamical behaviours of various classes.



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 30 of 33

8. Conclusion

With the help of new recently introduced concepts of fractional calculus, we have inves-
tigated a determinacy mathematical model of COVID-19. Because such modeling is very
effective if the involve parameters and functions are suitably defined. Additionally, if the
parameters can be accurately approximated, epidemic modeling is a very useful method
for estimating the state of the COVID-19 worldwide pandemic. Keeping these things in
mind, we have revisited a mathematical model [29] under the piecewise fractional order
derivative. Fundamental results including positivity, boundedness and equilibrium points
were deduced. Also reproductive numbers was calculated. Here, we have investigated sen-
sitivity analysis of the model via reproductive number analysis. We observe that increase
the value value of percentage upto some extent can increase or decrease the sensitivity
index by that percent. The existence theory is an important tool to be investigated for
dynamical problem. Here, we have used Schauder’s and Banach fixed point theorems to
establish a detailed analysis for existence theory of solution to the proposed model. Addi-
tional stability is important consequence and for dynamical systems two kinds of stability
be studied. One is devoted to stability of equilibrium points which is global or local in
nature. Therefore, we have attempted to deduce some results related to local asymptot-
ically stability of DFE point. It was deduced that DFE point is asymptotically stable if
R0 < 1. On the other hand Laypunov-Volterra theory was used to show some results for
global stability of endemic equilibrium point. The second stability results were devoted
to numerical type. Hence, we have used UH concepts to deduce such results. Finally, a
sophisticated numerical algorithm was derived by using the RK2 method. We have used
this scheme to simulate our results graphically. We see that crossover behaviors can seen
in each compartment after t1 in their graphical presentation.

In the future, we can investigate the mentioned models from different perspectives like:

• By using fractals fractional calculus the considered model has not studied.

• On applying the traditional tools of fractional calculus in terms of conformable can
be investigated.

• We can investigate the model by using non-singular type derivatives.

• Stochastic form of this model has not yet studied.

• In the future, if we include the harmonic mean type nonlinear incidence rate to
properly investigate the disease dynamics.

In the future, we will extend the numerical schemes introduced in [2–4] for piecewise
derivatives of fractional orders using different kinds of definitions.

Acknowledgements

We are thankful to Prince Sultan University for support through TAS research lab.



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 31 of 33

References

[1] R. P. Agarwal. Fixed Point Theory and Applications. Cambridge University Press,
2001.

[2] N. Almutairi and S. Saber. Chaos control and numerical solution of time-varying
fractional newton-leipnik system using fractional atangana-baleanu derivatives. AIMS
Mathematics, 8(11):25863–25887, 2023.

[3] N. Almutairi and S. Saber. On chaos control of nonlinear fractional newton-leipnik
system via fractional caputo-fabrizio derivatives. Scientific Reports, 13(1):22726,
2023.

[4] N. Almutairi and S. Saber. Application of a time-fractal fractional derivative with a
power-law kernel to the burke-shaw system based on newton’s interpolation polyno-
mials. MethodsX, 12:102510, 2024.

[5] M. H. Aqlan, A. Alsaedi, B. Ahmad, and J. J. Nieto. Mathematical modeling and
prediction in infectious disease epidemiology. Open Mathematics, 14(1):723–735, 2016.

[6] A. Atangana. Modelling the spread of covid-19 with new fractal-fractional operators:
Can the lockdown save mankind before vaccination? Chaos, Solitons Fractals,
136:109860, 2020.

[7] A. Atangana and S. I. Araz. New concept in calculus: Piecewise differential and
integral operators. Chaos, Solitons Fractals, 145:110638, 2021.

[8] A. Atangana and S. I. Araz. Rhythmic behaviors of the human heart with piecewise
derivative. Mathematical Biosciences and Engineering, 19:3091–3109, 2022.

[9] A. Atangana and S. I. Araz. Piecewise differential equations: Theory, methods and
applications. AIMS Mathematics, 8(7):15352–15382, 2023.

[10] N. Bacaer. A Short History of Mathematical Population Dynamics. Springer, London,
2011.

[11] F. Brauer, C. Castillo-Chavez, and Z. Feng. Mathematical Models in Epidemiology.
Springer, New York, 2019.

[12] C. Chen, M. Bohner, and B. Jia. Ulam-hyers stability of caputo fractional difference
equations. Mathematical Methods in the Applied Sciences, 42(18):7461–7470, 2019.

[13] A. K. Cheng. Mathematical modelling and real life problem solving. Association of
Mathematics Educators, 2009.

[14] Y. J. Cho. Survey on metric fixed point theory and applications. Advances in Real
and Complex Analysis with Applications, pages 183–241, 2017.

[15] M. Dalir and M. Bashour. Applications of fractional calculus. Applied Mathematical
Sciences, 4(21):1021–1032, 2011.

[16] L. Edelstein-Keshet. Mathematical Models in Biology. SIAM, 2004.
[17] A. Granas and J. Dugundji. Fixed Point Theory, volume 14. Springer, New York,

2003.
[18] A. B. Gumel, S. Ruan, T. Day, J. Watmough, F. Brauer, P. van den Driessche,

D. Gabrielson, C. Bowman, M. E. Alexander, S. Ardal, and J. Wu. Modelling strate-
gies for controlling sars outbreaks. Proceedings of the Royal Society of London. Series
B: Biological Sciences, 1554:2223–2232, 2004.



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 32 of 33

[19] A. Huppert and G. Katriel. Mathematical modeling and prediction in infectious
disease epidemiology. Clinical Microbiology and Infection, 19(11):999–1005, 2013.

[20] R. M. Jena. Piecewise concept in fractional models. Academic Press, New York, 2024.
[21] H. Khan, S. Ahmad, J. Alzabut, and A. T. Azar. A generalized coupled system of

fractional differential equations with application to finite time sliding mode control
for leukemia therapy. Chaos, Solitons Fractals, 174:113901, 2023.

[22] H. Khan, J. Alzabut, W. F. Alfwzan, and H. Gulzar. Nonlinear dynamics of a piece-
wise modified abc fractional-order leukemia model with symmetric numerical simula-
tions. Symmetry, 15(7):1338, 2023.

[23] H. Khan, J. Alzabut, J. F. Gomez-Aguilar, and R. P. Agarwal. Piecewise mabc
fractional derivative with an application. AIMS Mathematics, 8(10):24345–24366,
2023.

[24] S. Khan, Z. A. Khan, H. Alrabaiah, and S. Zeb. On using piecewise fractional differ-
ential operator to study a dynamical system. Axioms, 12(3):292, 2023.

[25] S. Khan, K. Shah, A. Debbouche, S. Zeb, and V. Antonov. Solvability and ulam-hyers
stability analysis for nonlinear piecewise fractional cancer dynamic systems. Physica
Scripta, 99(2):025225, 2024.

[26] V. Kiryakova. The special functions of fractional calculus as generalized fractional cal-
culus operators of some basic functions. Computers Mathematics with Applications,
59(3):1128–1141, 2020.

[27] R. L. Magin. Fractional calculus in bioengineering: A tool to model complex dynam-
ics. In Proceedings of the 13th International Carpathian Control Conference (ICCC),
pages 464–469. IEEE, 2012.

[28] M. Meerschaert. Mathematical Modeling. Academic Press, New York, 2013.
[29] C. Photphanloet, S. Ritraksa, S. E. Shuaib, A. Intarasit, and P. Riyapan. A com-

partmental model for assessing effects of covid-19 vaccination in thailand. Universal
Journal of Public Health, 10(6):596–605, 2022.

[30] S. S. Redhwan, M. Han, M. A. Almalahi, M. A. Alyami, M. Alsulami, and N. Al-
ghamdi. Piecewise implicit coupled system under abc fractional differential equations
with variable order. AIMS Mathematics, 9(6):15303–15324, 2024.

[31] M. Roccetti. Drawing a parallel between the trend of confirmed covid-19 deaths in
the winters of 2022/2023 and 2023/2024 in italy, with a prediction. Mathematical
Biosciences and Engineering, 21(3):3742–3754, 2024.

[32] B. Ross. A brief history and exposition of the fundamental theory of fractional
calculus. In Fractional Calculus and Its Applications. Springer, Berlin, 1975.

[33] I. A. Rus. Ulam stability of the operatorial equations. Springer, New York, 2011.
[34] M. A. Sanchez and S. M. Blower. Uncertainty and sensitivity analysis of the basic

reproductive rate: Tuberculosis as an example. American Journal of Epidemiology,
145(12):1127–1137, 1997.

[35] F. Schaffner and A. Mathis. Dengue and dengue vectors in the who european re-
gion: past, present, and scenarios for the future. The Lancet Infectious Diseases,
14(12):1271–1280, 2014.

[36] K. Shah, T. Abdeljawad, and A. Ali. Mathematical analysis of the cauchy type



S. Naz et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5861 33 of 33

dynamical system under piecewise equations with caputo fractional derivative. Chaos,
Solitons Fractals, 161:112356, 2022.

[37] K. Shah, A. Khan, B. Abdalla, T. Abdeljawad, and K. A. Khan. A mathematical
model for nipah virus disease by using piecewise fractional order caputo derivative.
Fractals, 32(02):2440013, 2024.

[38] K. Shah, A. Khan, B. Abdalla, T. Abdeljawad, and K. A. Khan. A mathematical
model for nipah virus disease by using piecewise fractional order caputo derivative.
Fractals, 32(2):2440013, 2024.

[39] K. Shah, P. Kumam, and I. Ullah. On ulam stability and multiplicity results to a
nonlinear coupled system with integral boundary conditions. Mathematics, 7(3):223,
2019.

[40] K. Shah, H. Naz, T. Abdeljawad, and B. Abdalla. Study of fractional order dynam-
ical system of viral infection disease under piecewise derivative. CMES-Computer
Modeling in Engineering Sciences, 136(1):921–941, 2023.

[41] N. H. Tuan, H. Mohammadi, and S. Rezapour. A mathematical model for covid-19
transmission by using the caputo fractional derivative. Chaos, Solitons Fractals,
140:110107, 2020.

[42] V. V. Uchaikin. Fractional Derivatives for Physicists and Engineers, volume 2.
Springer, Berlin, 2013.

[43] C. Vargas-De-León. Volterra-type lyapunov functions for fractional-order epidemic
systems. Communications in Nonlinear Science and Numerical Simulation, 24(1-
3):75–85, 2015.

[44] J. T. Wu, K. Leung, and G. M. Leung. Nowcasting and forecasting the potential
domestic and international spread of the 2019-ncov outbreak originating in wuhan,
china: a modelling study. The Lancet, 395(10225):689–697, 2020.

[45] A. Zeb, A. Atangana, Z. A. Khan, and S. Djillali. A robust study of a piece-
wise fractional order covid-19 mathematical model. Alexandria Engineering Journal,
61(7):5649–5665, 2022.

[46] Y. Zhang, X. Yu, H. Sun, G. R. Tick, W. Wei, and B. Jin. Applicability of time
fractional derivative models for simulating the dynamics and mitigation scenarios of
covid-19. Chaos, Solitons Fractals, 138:109959, 2020.

[47] S. Zhao, Q. Lin, J. Ran, S. S. Musa, G. Yang, W. Wang, Y. Lou, D. Gao, L. Yang,
D. He, and M. H. Wang. Preliminary estimation of the basic reproduction number
of novel coronavirus (2019-ncov) in china, from 2019 to 2020: A data-driven analysis
in the early phase of the outbreak. International Journal of Infectious Diseases,
92:214–217, 2020.

[48] Y. Zhou. Basic Theory of Fractional Differential Equations. World Scientific, Singa-
pore, 2010.