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

Fractional Discrete Modeling and Numerical Analysis of
Hepatitis B and C Dynamics with Neural Network

Approximation

Sajid Ali Shah1, Nasir Rehman2, Sadique Ahmad3,∗,Mohammed A.Elaffendi3,
Abdelhamied Ashraf Ateya3

1 Department of Zoology Government Degree College Thana Malakand, Khyber Pakhtunkhwa,
Pakistan
2 Department of Mathematics, Allama Iqbal Open University Islamabad, Pakistan
3 EIAS Data Science and BlockChain Laboratory, College of Computer and Information
Sciences, Prince Sultan University, Riyadh 11586, Saudi Arabia

Abstract. This paper is devoted to establish a details analysis for the dynamics of hepatitis B
and C through some mathematical difference equations. Serious viral illnesses that mainly affect
the liver, hepatitis B and C can cause inflammation and even damage to the liver. Hepatitis B can
spread by bodily fluids such as semen or vaginal fluid, even though both are blood-borne. While
hepatitis C can be cured with medication, hepatitis B can be prevented with vaccine. We formu-
late a mathematical model of the mentioned disease through fractional difference equations using
the Caputo fractional difference operator (CFDO). Existence theory, Ulam-Hyers (UH) stability
and numerical investigation are deduced. Further, analysis based on artificial intelligence (AI) is
considered. We have used neural networks (NNs) based investigations to conduct our study. Var-
ious graphical illustrations are given to demonstrate our results. For our analysis we use Matlab
software with 2023 version.

2020 Mathematics Subject Classifications: 26A33, 34A08, 68T07, 92B20

Key Words and Phrases: Viral infections, Difference equations, Numerical investigation, Neural
networks.

1. Introduction

A well-known area of study in the biological sciences, epidemiology looks at health,
disease, and associated variables at the population level. The science of epidemiology
is primarily concerned with human populations, as evidenced by the fact that the term

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i3.6481

Email addresses: shahsajidali442@gmail.com (S. A. Shah),
nasir.rehman@aiou.edu.pk (N. Rehman),saahmad@psu.edu.sa (S. Ahmad),
Affendi@psu.edu.sa (M. A. Affendi), aateya@psu.edu.sa (A. A.ateya)

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



S. A.Shah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6481 2 of 19

”epidemiology” is formed from three Greek words: epi (upon), demos (public), and logos
(research) [1]. Nowadays, there is a wealth of opportunities to use mathematical research
in the biological sciences to the study of infectious diseases. Both human and non-human
populations’ rates of illness and mortality are greatly influenced by infectious diseases.
As a result, both mathematicians and biologists are interested in the potential of formal
models to alleviate these pressures [2]. In the field of mathematical biology, numerous
scholars have devoted their time to investigating the transmission of infectious illnesses
using a variety of well-established mathematical models [3]. Throughout history, infectious
diseases have had a significant impact on people’s lives. In many parts of the world,
infectious diseases that manifest as outbreaks cause pandemic that claim millions of lives.
For example, the 1918 influenza virus [4], commonly referred to as ”the Spanish Flu,”
produced the greatest pandemic in recorded history before COVID-19 [5]. Over 50 million
people died globally during the aforementioned outbreak. Since the spread of susceptible
and ill individuals might infect more people and potentially result in their death, infectious
diseases have long posed a threat to human safety [6].

Globally, between 64 and 103 million people have a chronic infection. Unsterile medical
procedures, iatrogenic infections, and inappropriate injectable drug use are key risk factors
for HCV infection in countries where it is extremely prevalent. In regions with high HCV
incidence, unsterile medical procedures (iatrogenic infections) and improper injectable
drug use are major risk factors for these blood-borne viral infections [7]. Hepatitis B
is a viral disease that can cause liver failure, cirrhosis, and liver cancer. More than
250 million persons worldwide are estimated to have a chronic Hepatitis B infection by
the World Health Organization (WHO) [8]. Hepatitis C is a virus that can cause liver
cancer, cirrhosis, and liver failure. Over 70 million people worldwide suffer from a chronic
Hepatitis C infection, according to the WHO [9]. Hepatitis B and C are major public
health concerns in Pakistan. According to a nationwide survey conducted by the Pakistan
Medical Research Council (PMRC) in 2017-18, the prevalence of Hepatitis B and C in the
general population aged 15 and older was estimated to be 2.5% and 4.8%, respectively
[10].

It is crucial to remember that when real-world issues are represented in terms of math-
ematical equations, they become more effective instruments for comprehending phenom-
ena. The study of mathematical epidemiology has drawn much more attention in the
last 200 years [11]. Numerous infectious diseases have been studied using the concept of
mathematical biology. Much research effort has been focused on the aforementioned field.
Kermack-McKendrick created a simple mathematical model called susceptible, infected,
and recovered (SIR). Only infection and removal events are included in this model, which
is adequate to explain a basic epidemic, including the threshold condition required for an
epidemic to start (see [12]). Mathematical modeling is a very clever tool that scientists
and engineers use to solve real-world problems [13]. Recently, fractional calculus has been
used in modeling real world problems interms of fractional difference equations (FDEs).
The literature has put forth a number of definitions of fractional-order operators, each
with unique advantages and disadvantages. As a logical progression of the integer-order
difference equation, one such operator is called a fractional-order discrete difference opera-



S. A.Shah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6481 3 of 19

tor. Since real-world data is usually in discrete form, these operators work incredibly well
for researching real-world phenomena since they produce output in the discrete form. Im-
age processing[14] , game theory[15], ecology[16], neural networks[17], epidemiology[18],
economics[19], and electrical engineering[20] are just a few of its many uses. Recently,
researchers have explored the area further and produced significant results, we refer to
[21, 22].

Because it makes it easier to evaluate experimental results and comprehend the bi-
ological mechanisms driving epidemic spread, mathematical modeling is a valuable tool
for studying virus dynamics. Theoretical biologists can better understand the processes
that contribute to the disease by using mathematical models to study the dynamics of the
Hepatitis B virus. Numerous models have been proposed to comprehend the dynamics of
Hepatitis B and C. Authors have studied the model [23] for Hepatitis B and C given by

Ḣ(t) = Λ− δH − ηHV,

İ(t) = ηHV − γI,

V̇ (t) = λI − µV,

H(0) = H0, I(0) = I0, V (0) = V0,

(1)

where H stands for healthy cells, I for infected and V for virus free cells density respec-
tively. Λ is the birth rate, η is efficacy rate of the process to infection, δ, µ and γ stand for
death rates respectively of healthy, infected and virus free cells respectively. In addition, λ
is the free virus cell production rate. We extend the model (1) under the Caputo fractional
difference operator as described by

C∆ξH(t) = Λ− δH(t+ ξ − 1)− ηH(t+ ξ − 1)V (t+ ξ − 1),
C∆ξI(t) = ηH(t+ ξ − 1)V (t+ ξ − 1)− γI(t+ ξ − 1),
C∆ξR(t) = λI(t+ ξ − 1)− µV (t+ ξ − 1),

(2)

where t ∈ [0, b]N0 . We describe the formulation (2) in schematic diagram as given in figure
1.



S. A.Shah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6481 4 of 19

Figure 1: Schematic diagram of model (2).

We use fixed point theorem to deduce the existence theory for model (2). Also, numer-
ical analysis against various fractional orders has been given. UH stability is also deduced
for the considered model. In addition, in recent time AI tools based NNs have attracted
attention to use in investigation of different epidemic models. Large amounts of data are
utilized to train data-driven neural networks, which can be used to build mathematical
models and solve inverse problems [24]. Data-driven neural networks can assist in improv-
ing the current models and equations and obtaining an equation that describes the process
being studied. Deep learning techniques can be used to do this. For instance, authors
[25] was suggested to use deep learning to forecast how many people in various parts of
the US would contract COVID-19 and pass away over the course of one to four weeks.
AI models based on deep neural networks and epidemic compartmental models based on
differential equations are effective tools for assessing and preventing the spread of viral
infectious disease like HBV and HCV. However, the difficulties in estimating parameters
restrict the potential of compartmental models, and AI models are unable to identify the
evolutionary pattern of said disease and are not explainable [26, 27]. The concerned NNs
are consisted on multi layers including input output layers. The structure of such NNs,
we can describe as in figure 3.



S. A.Shah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6481 5 of 19

Figure 2: Structure of NNs.

With the help of mentioned NNs, we deduce some analysis for the concerned model
by computing the root mean square (RMS) and mean square error (MSE) as well as re-
gression coefficient R. Neural networks learn, optimize, and forecast using a variety of
algorithms. Back propagation, gradient descent, and its variations, such as Adam and
Stochastic Gradient Descent (SGD), are important algorithms. Convolutional neural net-
works (CNNs), recurrent neural networks (RNNs), and long short-term memory (LSTM)
networks are further significant methods. Neural models are trained using the Liveners-
Marquardt backward propagation algorithm (LMBA). Recently, researchers have used
NNs and fractional calculus to investigate various problems, we refer to [28–30].

2. Basic Results

Some basic results are demonstrated below.

Definition 1. [18] For function χ(t), the ξ-CFDO is defined by

C∆ξ
aχ(t) = ∆−(n−ξ)

a ∆nχ(t) =
1

Γ(n− ξ)

t−(n−ξ)∑
s=a

(t− s− 1)n−ξ−1∆nχ(s), (3)

where t ∈ Na+n−ξ,n = ⌈ξ⌉+1 and ξ > 0 with ξ /∈ N and (t− 1− s)(n−ξ−1), while the n-the
integer difference operator denoted by ∆nχ(s) which can be expressed as follows:

(t− s− 1)(n−ξ−1) =
Γ(t− s)

Γ(t− s− n + ξ + 1)
, (4)

and

∆nχ(t) = ∆(∆n−1χ(t)) =

n∑
k=0

(
n

k

)
(−1)n−kχ(t+ k), t ∈ Na. (5)



S. A.Shah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6481 6 of 19

Remark 1. For n = 1, the ξ-CFDO is defined as:

C∆ξ
af(t) = ∆−(1−ξ)

a ∆χ(t) =
1

Γ(1− ξ)

t−(1−ξ)∑
s=a

(t− s− 1)−ξ∆χ(s), (6)

where t ∈ Na+1−ξ.

Definition 2. [18] With the mapping χ : Na → R, let χ(t) be a function defined on Na.
The following is the definition of the ξ-the fractional sum of f in this context:

∆−ξχ(t) =
1

Γ(ξ)

t−ξ∑
s=a

(t− s− 1)ξ−1χ(s), (7)

where ξ > 0, t ∈ Nb+ξ = {b+ ξ, b+ ξ + 1, . . . } is a time scale and b ∈ R.

Lemma 1. [19] Let ξ > 0 and χ be defined on Na, then:

∆−ξC∆ξχ(t) = χ(t)−
∑n−1

k=0
(t−b)(k)

k! ∆kh(b)

= f(t) + c0 + c1t+ · · ·+ cnt
n−1,

such that n satisfies n ≥ ξ and ci ∈ R, i = 1, 2, . . . ,n− 1.

Lemma 2. [19] The following equality holds:

t−ξ∑
s=0

(t− s− 1)(ξ−1) =
Γ(t+ 1)

ξΓ(t− ξ + 1)
. (8)

The proposed model (2) can be written as follows:

H(t) = H(0) +
1

Γ(ξ)

t−ξ∑
s=1−ξ

(t− s− 1)ξ−1(Λ− ηH(t+ ξ − 1)V (t+ ξ − 1)− δH(t+ ξ − 1))

I(t) = I(0) +
1

Γ(ξ)

t−ξ∑
s=1−ξ

(t− s− 1)ξ−1(ηH(t+ ξ − 1)V (t+ ξ − 1)− γI(t+ ξ − 1)),

V (t) = V (0) +
1

Γ(ξ)

t−ξ∑
s=1−ξ

(t− s− 1)ξ−1(λI(t+ ξ − 1)− µV (t+ ξ − 1).

(9)

3. Existence and Uniqueness

Take model (9) in the following form
χ1(t+ ξ − 1, H, I, V ) =Λ− ηH(t+ ξ − 1)V (t+ ξ − 1)− δH(t+ ξ − 1),

χ2(t+ ξ − 1, H, I, V ) =ηH(t+ ξ − 1)V (t+ ξ − 1)− γI(t+ ξ − 1),

χ3(t+ ξ − 1, H, I, V ) = λI(t+ ξ − 1)− µV (t+ ξ − 1),

(10)



S. A.Shah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6481 7 of 19

then equation (9) may be described as

g(t) = g(0) +
1

Γ(ξ)

t−ξ∑
s=0

(t− s− 1)ξ−1F (s + ξ − 1, g(s + ξ − 1)), (11)

while

g(t) =


H(t)

I(t)

V (t)

, g(0) =


H0

I0

V0

, F (t+ ξ− 1, g(t+ ξ− 1)) =


χ1(t+ ξ − 1, H, I, V )

χ2(t+ ξ − 1, H, I, V )

χ3(t+ ξ − 1, H, I, V )

(12)

The following fixed-point theorems must be applied in order to prove the existence and
uniqueness of solutions for the system (9).

Theorem 1. [27] A contraction mapping established on a complete metric space has ex-
actly one fixed point, according to the unique fixed point feature.

Theorem 2. [27] . This means that ψ has at least one fixed point in D if ψ : D ⊆ Rn →
D ⊆ Rn is a continuous mapping defined on a nonempty, bounded, closed, and convex set
D.

We define the Banach space by X = X1 × X2 × X3 such that Xi, (i = 1, 2, 3) are

sequences of real numbers (H, I, V ) = {H(t), I(t), V (t)}b+ξt=ξ−1 under the norm ∥g∥ =
∥(H, I, V )∥ = supt∈[ξ−1,b+ξ]Nξ−1

[|H(t)|+ |I(t)|+ |V (t)|]. Then ψX → X such that

(ψg)(t) = g(0) +
1

Γ(ξ)

t−ξ∑
s=0

(t− s− 1)ξ−1F (s + ξ − 1, g(s + ξ − 1)). (13)

Define Dψ = max
{

Γ(t+1)
Γ(ξ+1)Γ(t−ξ+1)

}
= Γ(b+ξ+1)

Γ(ξ+1)Γ(b+1) .

Theorem 3. Let

(C1) : We have bounded function F : [ξ − 1, b + ξ]Nξ−1
→ R with

∣∣F (t, g)∣∣ ≤ F(t)
∣∣g∣∣ at

every g ∈ X.

The system (2) has at least one solution at X, if F∗ < Dψ, with F∗ = max{F(t) : t ∈
[ξ − 1, b+ ξ]Nξ−1

}.

Proof. Take P > 0 such that U =
{
g
∣∣[ξ − 1, b + ξ]Nξ−1

→ R, ∥g∥ ≤ M
}
, we need to

show that ψ maps U into U . For any g ∈ U , we have

∣∣(ψg)(t)∣∣ ≤ F(t)

Γ(ξ)

t−ξ∑
s=0

(t− s− 1)ξ−1
∣∣g(s + ξ − 1)

∣∣,
≤ F(t)∥g∥

Γ(ξ)

t−ξ∑
s=0

(t− s− 1)ξ−1,



S. A.Shah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6481 8 of 19

≤ DψF∗P.

Thus |(ψg)(t)| < P, then ∥(ψg)∥ ≤ P implies that ψ maps U to itself. On using Theorem
(2), ψ has at least one solution. Hence system (2) has also at least one solution.

Theorem 4. Let

(C2) : If there exists constant K > 0 such that
∣∣F (t, g) − F (t, ḡ)

∣∣ ≤ K
∣∣g − ḡ

∣∣ at all
t ∈ [ξ − 1, b+ ξ]Nξ−1

and g, ḡ ∈ X.

The system (2) has a unique solution on X if K < Dψ.

Proof. If g, ḡ ∈ X, then for every t ∈ [ξ − 1, b+ ξ]Nξ−1
such that

∣∣(ψg)(t)− (ψḡ)(t)
∣∣ ≤ 1

Γ(ξ)

t−ξ∑
s=0

(t− s− 1)ξ−1
∣∣F (s + ξ − 1, g(s + ξ − 1))

−G(s + ξ − 1, g(s + ξ − 1))
∣∣,

≤ F
Γ(ξ)

t−ξ∑
s=0

(t− s− 1)ξ−1
∣∣g(s + ξ − 1))− ḡ(s + ξ − 1))

∣∣,
≤ F∥g − ḡ∥

Γ(ξ)

t−ξ∑
s=0

(t− s− 1)ξ−1,

≤ DψF∥g − ḡ∥.

Thus ∥(ψg) − (ψḡ)∥ ≤ ∥g − ḡ∥, Since ψ is a contraction, there is a unique fixed point
of ψ that is also a unique solution of the system (2) according to the Banach fixed-point
theorem (1).

Here, we deduce some results for stability analysis.

|C∆ξg(t)− F (s + ξ − 1, g(s + ξ − 1))| ≤ ϵ, t ∈ [0, b]N0 . (14)

Definition 3. [28] The problem (2) is UH stable for any ϵ > 0, if we have a constant
Dg > 0 which obeys the inequality (14) such that

|g(t)− ḡ(t)| ≤ Dgϵ, t ∈ [ξ − 1, b+ ξ]Nξ−1
. (15)

The generalization of equation (2) If the constant Dgϵ in inequality (14) is replaced by
the function ϕ(ϵ), UH stability is maintained, where ϕ ∈ C(R+,R+) and ϕ(0) = 0.

Remark 2. If g ∈ X is a solution to inequality (14) and there exists Θ : [ξ−1, b+ξ]Nξ−1
→

R, ∋

(i) |Θ(s + ξ − 1)| ≤ ϵ, t ∈ [0, b]N0 ,

(ii) C∆ξg(t) = F (s + ξ − 1, g(s + ξ − 1)) + Θ(s + ξ − 1), t ∈ [0, b]N0 .



S. A.Shah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6481 9 of 19

Theorem 5. In view of (C2), the suggested model (2) is UH stable if K < D−1
ψ .

Proof. Take g, ḡ ∈ X be the solution of equation (2) using the Remark 2, one has

g(t) = g(0) +
1

Γ(ξ)

t−ξ∑
s=0

(t− s− 1)ξ−1F (s + ξ − 1, g(s + ξ − 1)) + Θ(s + ξ − 1). (16)

In view of Remark 2, for t ∈ [ξ − 1, b+ ξ]Nξ−1
}, we get∥∥∥g(t)− (

g(0) +
1

Γ(ξ)

t−ξ∑
s=0

(t− s− 1)ξ−1F (s + ξ − 1, g(s + ξ − 1))
)∥∥∥ ≤ Dψϵ. (17)

For t ∈ [ξ − 1, b+ ξ]Nξ−1
, one has

|g(t)− ḡ(t)| =
∣∣∣g(t)− (

g(0) +
1

Γ(ξ)

t−ξ∑
s=0

(t− s− 1)ξ−1F (s + ξ − 1, g(s + ξ − 1))
)∣∣∣,

≤ Dψϵ+
1

Γ(ξ)

t−ξ∑
s=0

(t− s− 1)ξ−1

× |F (s + ξ − 1, g(s + ξ − 1))− F (s + ξ − 1, ḡ(s + ξ − 1))|,

≤ Dψϵ+
1

Γ(ξ)

t−ξ∑
s=0

(t− s− 1)ξ−1K∥g − ḡ∥ (18)

which further means

∥g − ḡ∥ ≤
[ Dψ

1−KDψ

]
ϵ = Dgϵ. (19)

Also set Dg =
Dψ

1−KDψ > 0, we conclude that solution of (2) is UH stable. It is possible

to determine that the solution of model (2) satisfies the requirements for generalized UH
stability by choosing ϕ(ϵ) = Dgϵ, where ϕ(0) = 0.

4. Qualitative Results of Model (2)

Here, we deduce some results related to positivity, equilibrium points and reproductive
number. Also, stability result is deduced. We see that from Eq.(9)

H(t) = H(0) +
1

Γ(ξ)

t−ξ∑
s=1−ξ

(t− s− 1)ξ−1(Λ− ηH(t+ ξ − 1)V (t+ ξ − 1)− δH(t+ ξ − 1)) > 0

I(t) = I(0) +
1

Γ(ξ)

t−ξ∑
s=1−ξ

(t− s− 1)ξ−1(ηH(t+ ξ − 1)V (t+ ξ − 1)− γI(t+ ξ − 1)) > 0,

V (t) = V (0) +
1

Γ(ξ)

t−ξ∑
s=1−ξ

(t− s− 1)ξ−1(λI(t+ ξ − 1)− µV (t+ ξ − 1) > 0,



S. A.Shah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6481 10 of 19

which implies that H(t), I(t), V (t) are positive for t > 0. The disease free E0 and
endemic E∗ equilibrium points are computed as follows: E0 = (H0, I0, V0) =

(
Λ
δ , 0, 0

)
,

and E∗ = (H∗, I∗, V ∗), where

H∗ =
γµ

ηλ
, I∗ =

µ(ληΛ− δγµ)

λ3η
, V ∗ =

ληΛ− δγµ

λ2η
.

Further, the reproductive number is computed as

R0 =
ηΛλ

δγµ
.

Theorem 6. The point E0 is asymptotically stable if R0 < 1.

Proof. From the Jacobian matrix J at E0, we have

det(J − ρI) =

∣∣∣∣∣∣
ρ+ δ 0 ηΛ

δ

0 ρ+ γ −ηΛ
δ

0 λ ρ+ µ

∣∣∣∣∣∣ = 0

which yields that ρ = −δ < 0, and ρ2 + a1ρ + a2 = 0, where a1 = γ + µ, a2 = γδµ−Λλη
η ,

clearly a1 > 0 and a2 > 0 if and only if R0 < 1.

Theorem 7. The point E∗ is locally asymptotically stable if R0 ≥ 1.

Proof. The endemic equilibrium point is given by E∗ = (1.0714285714, 0.0849560567, 0.0004757539).
The characteristic equations becomes det(J(E∗) − ρI) = 0. Upon computation, we get
eigens values as

ρ1 = −0.0345678−0.0543217ι, ρ2 = −0.0453276+0.079805412ι, ρ3 = −0.005123984+0ι.

We see that all eigen’s values have negative real parts. Thus the said equilibrium point is
locally asymptotically stable.

5. Numerical Analysis

We reset the system (9) at s + ξ = k to get the following system

H(n) = H(0) +
1

Γ(ξ)

n∑
k=1

Γ(n− k + ξ)

Γ(n− k + 1)
(Λ− ηH(k − 1)I(k − 1)− δH(k − 1)),

I(n) = I(0) +
1

Γ(ξ)

n∑
k=1

Γ(n− k + ξ)

Γ(n− k + 1)
(ηH(k − 1)I(k − 1)− γI(k − 1)),

V (n) = V (0) +
1

Γ(ξ)

n∑
k=1

Γ(n− k + ϕ)

Γ(n− k + 1)
(λI(k − 1)− µV (k − 1)).

(20)



S. A.Shah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6481 11 of 19

We take the initial values in millions as H(0) = 10, I(0) = 2, V (0) = 5, and λ =
0.0003,Λ = 0.02, δ = 0.0001, η = 0.0025, γ = 0.05, µ = 0.0056. We present the nu-
merical interpretation for the three classes in figures 3 and figure 4 for different fractional
order values.

0 50 100 150 200 250 300 350 400 450

time t (Days)

1

2

3

4

5

6

7

8

9

H
ea

lt
h

y 
ce

lls

(a)

0.75
0.85
0.95
1.0

0 50 100 150 200 250 300 350 400 450 500

time t (Days)

0

0.5

1

1.5

2

2.5

In
fe

ct
ed

 c
el

ls

(b)

0.75
0.85
0.95
1.0

0 50 100 150 200 250 300 350 400 450 500

time t (Days)

5

5.2

5.4

5.6

5.8

6

6.2

6.4

V
ir

u
s 

fr
ee

 c
el

ls

(c)

0.75
0.85
0.95
1.0

Figure 3: Graphical presentations of approximate solutions using different fractional orders
in [0.75, 1.0] for model (2) .



S. A.Shah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6481 12 of 19

0 50 100 150 200 250 300 350 400 450 500

time t (Days)

1

2

3

4

5

6

7

8

9

10

H
ea

lt
h

y 
ce

lls

(a)

0.35
0.45
0.55
0.65

0 50 100 150 200 250 300 350 400 450 500

time t (Days)

0

0.5

1

1.5

2

2.5

In
fe

ct
ed

 c
el

ls

(b)

0.35
0.45
0.55
0.65

0 50 100 150 200 250 300 350 400 450 500

time t (Days)

5

5.5

6

6.5

V
ir

u
s 

fr
ee

 c
el

ls

(c)

0.35
0.45
0.55
0.65

Figure 4: Graphical presentations of approximate solutions using different fractional orders
in [0.35, 0.65] for model (2) .

In above both figures 3, 4, we see that density of healthy cells is decreasing with
different fractional order until achieve stable position. In the same way, the infected cells
density is first rising exponentially and then it goes on decreasing with the passage of
time under various fractional order. The density of virus free cells is gradually rising until
reaches to a stable position by using different fractional order.

6. NNs Analysis of Model (2)

Here, we demonstrate the analysis based on AI NNs tools. We compute the MSE,
RMSE and regression coefficient. The circuit diagram of the concerned NNs is given in
figure 5.

Figure 5: Diagram of neural network that we have used.



S. A.Shah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6481 13 of 19

Onward, we presents the graphical illustrations of various data including all data, test
data, validation and train data. Also the function fit for the given compartments have been
presented. In figure 6 sub plots (a)− (d), we have presented all data, test data, validation
of data and train data for the density of healthy cells using adopted NNs. Further, in sub
plot (e) of figure 6, we have presented the absolute error plots between the numerical and
and NNs data analysis, we see that error is very low which shows the best performance
of our methodology. Also, the function fits corresponding to various data has given in
subplot (f) of figure 6. Further, we have used NNs with two hidden layers of 20 neurons.
For other units, we give table 1.

0 20 40 60 80 100 120
6.5

7

7.5

8

8.5

9

9.5

10
All Data (H)

Outputs
Targets

6.5 7 7.5 8 8.5 9 9.5 10
6.5

7

7.5

8

8.5

9

9.5

10
R = 1.0

0 20 40 60 80 100 120
-0.06

-0.04

-0.02

0

0.02

0.04
MSE = 6.7725e-05, RMSE = 0.0082295

Error

 = -0.00092669,  = 0.008218

-0.05 -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03
0

5

10

15

20

25

30

(a)

0 10 20 30 40 50 60 70 80
6.5

7

7.5

8

8.5

9

9.5

10
Train Data (H)

Outputs
Targets

6.5 7 7.5 8 8.5 9 9.5 10
6.5

7

7.5

8

8.5

9

9.5

10
R = 1

0 10 20 30 40 50 60 70 80
-6

-4

-2

0

2

4

6

8
10-3 MSE = 2.5751e-06, RMSE = 0.0016047

Error

 = -0.00042583,  = 0.0015582

-4 -2 0 2 4 6

10-3

0

2

4

6

8

10

(b)

0 5 10 15
6.5

7

7.5

8

8.5

9

9.5

10
Validation Data (H)

Outputs
Targets

6.5 7 7.5 8 8.5 9 9.5 10
6.5

7

7.5

8

8.5

9

9.5

10
R = 1.0

0 5 10 15
-0.06

-0.05

-0.04

-0.03

-0.02

-0.01

0

0.01

0.02
MSE = 0.00030788, RMSE = 0.017547

Error

 = -0.0029459,  = 0.017905

-0.06 -0.04 -0.02 0 0.02 0.04 0.06
0

0.5

1

1.5

2

2.5

3

(c)

0 5 10 15
6.5

7

7.5

8

8.5

9

9.5

10
Test Data (H)

Outputs
Targets

6.5 7 7.5 8 8.5 9 9.5 10
6.5

7

7.5

8

8.5

9

9.5

10
R = 1.0

0 5 10 15
-0.02

-0.01

0

0.01

0.02

0.03
MSE = 0.00013595, RMSE = 0.01166

Error

 = -0.0012782,  = 0.011996

-0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04
0

0.5

1

1.5

2

(d)

0 5 10 15 20 25 30 35 40 45 50
t

10-7

10-6

10-5

10-4

10-3

10-2

10-1

A
b

so
lu

te
 e

rr
o

rs
 H

(t
)

0.75
0.85
0.95
1.00

(e)

0 5 10 15 20 25 30 35 40 45 50
6.5

7

7.5

8

8.5

9

9.5

10

O
u

tp
u

t 
an

d
 T

ar
g

et

Function Fit for Output Element H

Training Targets
Training Outputs
Validation Targets
Validation Outputs
Test Targets
Test Outputs
Errors
Fit

Input
-0.06

-0.04

-0.02

0

0.02

0.04

E
rr

o
r

Targets - Outputs

(f)

Figure 6: Graphical presentation of various data results for class H using NNs (a) all data,
(b) train data, (c) test data (d) validation (e) absolute error using NNs and (f) function
fit performance.



S. A.Shah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6481 14 of 19

Unite Initial value stopped value target value

Epoch 0 11 1000
Elapsed time - 00:00:00 -
Performance 24.9 4.06e− 11 0
Gradient 41.7 6.53e− 05 1e− 07
standard deviation 0.001 1e− 11 1e+ 11
Validation checks 0 6 6

Table 1: Training progress presentation for class H.

In figure 7 sub plots (a)− (d), we have presented all data, test data, validation of data
and train data for the density of infected cells using adopted NNs. Further, in sub plot
(e) of figure 7, we have presented the absolute error plots between the numerical and and
NNs data analysis, we see that error is very low which shows the best performance of our
methodology. Also, the function fit of class I corresponding to various data has given in
subplot (f) of figure 7. The NNs used contain two hidden layers with 20 neurons, while
other units are presented in table 2.



S. A.Shah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6481 15 of 19

0 20 40 60 80 100 120
2

2.05

2.1

2.15

2.2
All Data (I)

Outputs
Targets

2 2.02 2.04 2.06 2.08 2.1 2.12 2.14 2.16 2.18
2

2.05

2.1

2.15

2.2
R = 1.0

0 20 40 60 80 100 120
-2

-1

0

1

2

3

4
10-3 MSE = 4.2709e-07, RMSE = 0.00065352

Error

 = 0.00014152,  = 0.00064119

-2 -1 0 1 2 3

10-3

0

10

20

30

40

50

60

70

(a)

0 10 20 30 40 50 60 70 80
2

2.05

2.1

2.15

2.2
Train Data (H)

Outputs
Targets

2 2.02 2.04 2.06 2.08 2.1 2.12 2.14 2.16 2.18
2

2.05

2.1

2.15

2.2
R = 1

0 10 20 30 40 50 60 70 80
-3

-2

-1

0

1

2

3
10-5 MSE = 5.9067e-11, RMSE = 7.6855e-06

Error

 = -2.0324e-10,  = 7.7402e-06

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

10-5

0

5

10

15

20

(b)

0 5 10 15
2.02

2.04

2.06

2.08

2.1

2.12

2.14

2.16

2.18
Validation Data (H)

Outputs
Targets

2.02 2.04 2.06 2.08 2.1 2.12 2.14 2.16 2.18
2.02

2.04

2.06

2.08

2.1

2.12

2.14

2.16

2.18
R = 1.0

0 5 10 15
-2

-1

0

1

2

3

4
10-3 MSE = 1.9669e-06, RMSE = 0.0014025

Error

 = 0.00036993,  = 0.0014003

-4 -3 -2 -1 0 1 2 3 4 5

10-3

0

0.5

1

1.5

2

2.5

3

(c)

0 5 10 15
2.06

2.08

2.1

2.12

2.14

2.16

2.18
Test Data (I)

Outputs
Targets

2.06 2.08 2.1 2.12 2.14 2.16 2.18
2.06

2.08

2.1

2.12

2.14

2.16

2.18
R = 1.0

0 5 10 15
-1

-0.5

0

0.5

1

1.5

2
10-3 MSE = 9.0856e-07, RMSE = 0.00095318

Error

 = 0.000583,  = 0.00078057

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3

10-3

0

0.5

1

1.5

2

(d)

0 5 10 15 20 25 30 35 40 45 50
t

10-9

10-8

10-7

10-6

10-5

10-4

10-3

10-2

A
b

so
lu

te
 e

rr
o

rs
 I(

t)

0.75
0.85
0.95
1.00

(e)

0 5 10 15 20 25 30 35 40 45 50
2

2.02

2.04

2.06

2.08

2.1

2.12

2.14

2.16

2.18

O
u

tp
u

t 
an

d
 T

ar
g

et

Function Fit for Output Element I

Training Targets
Training Outputs
Validation Targets
Validation Outputs
Test Targets
Test Outputs
Errors
Fit

Input
-2

0

2

4

E
rr

o
r

10-3

Targets - Outputs

(f)

Figure 7: Graphical presentation of various data results for class H using NNs (a) all data,
(b) train data, (c) test data (d) validation (e) absolute error using NNs and (f) function
fit performance.

Unite Initial value stopped value target value

Epoch 0 17 1000
Elapsed time - 00:00:00 -
Performance 33.7 2.96e− 09 0
Gradient 56 9.73e− 08 1e− 07
standard deviation 0.001 1e− 09 1e+ 10
Validation checks 0 0 6

Table 2: Training progress presentation for class I.



S. A.Shah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6481 16 of 19

In figure 8 sub plots (a)− (d), we have presented all data, test data, validation of data
and train data for the density of virus free cells using adopted NNs. Further, in sub plot
(e) of figure 8, we have presented the absolute error plots between the numerical and and
NNs data analysis, we see that error is very low which shows the best performance of our
methodology. Also, the function fit for virus free cells density corresponding to various
data has given in subplot (f) of figure 8.

0 20 40 60 80 100 120
5

5.05

5.1

5.15

5.2

5.25

5.3

5.35

5.4
All Data (V)

Outputs
Targets

5 5.05 5.1 5.15 5.2 5.25 5.3 5.35 5.4
5

5.05

5.1

5.15

5.2

5.25

5.3

5.35

5.4
R = 1.0

0 20 40 60 80 100 120
-2

-1

0

1

2

3

4
10-3 MSE = 6.1351e-07, RMSE = 0.00078327

Error

 = 0.00013831,  = 0.0007748

-3 -2 -1 0 1 2 3

10-3

0

10

20

30

40

50

60

70

80

(a)

0 10 20 30 40 50 60 70 80
5

5.05

5.1

5.15

5.2

5.25

5.3

5.35

5.4
Train Data (V)

Outputs
Targets

5 5.05 5.1 5.15 5.2 5.25 5.3 5.35 5.4
5

5.05

5.1

5.15

5.2

5.25

5.3

5.35

5.4
R = 1

0 10 20 30 40 50 60 70 80
-1

-0.5

0

0.5

1

1.5
10-6 MSE = 2.6124e-13, RMSE = 5.1112e-07

Error

 = 1.4706e-08,  = 5.1455e-07

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

10-6

0

0.5

1

1.5

2

2.5

3

3.5

4

(b)

0 5 10 15
5

5.05

5.1

5.15

5.2

5.25

5.3

5.35
Validation Data (V)

Outputs
Targets

5 5.05 5.1 5.15 5.2 5.25 5.3 5.35
5

5.05

5.1

5.15

5.2

5.25

5.3

5.35
R =1.0

0 5 10 15
-2

-1

0

1

2

3
10-3 MSE = 1.9495e-06, RMSE = 0.0013963

Error

 = 0.00040112,  = 0.0013843

-4 -3 -2 -1 0 1 2 3 4 5

10-3

0

0.5

1

1.5

2

(c)

0 5 10 15
5

5.1

5.2

5.3

5.4
Test Data (H)

Outputs
Targets

5 5.05 5.1 5.15 5.2 5.25 5.3 5.35 5.4
5

5.1

5.2

5.3

5.4
R = 1.0

0 5 10 15
-2

-1

0

1

2

3

4
10-3 MSE = 2.1814e-06, RMSE = 0.001477

Error

 = 0.00053007,  = 0.001427

-4 -3 -2 -1 0 1 2 3 4 5

10-3

0

0.5

1

1.5

2

(d)

0 5 10 15 20 25 30 35 40 45 50
t

10-9

10-8

10-7

10-6

10-5

10-4

10-3

10-2

A
b

so
lu

te
 e

rr
o

rs
 V

(t
)

0.75
0.85
0.95
1.00

(e)

0 5 10 15 20 25 30 35 40 45 50
5

5.05

5.1

5.15

5.2

5.25

5.3

5.35

5.4

O
u

tp
u

t 
an

d
 T

ar
g

et

Function Fit for Output Element V

Training Targets
Training Outputs
Validation Targets
Validation Outputs
Test Targets
Test Outputs
Errors
Fit

Input
-2

0

2

4

E
rr

o
r

10-3

Targets - Outputs

(f)

Figure 8: Graphical presentation of various data results for class V using NNs (a) all data,
(b) train data, (c) test data (d) validation (e) absolute error using NNs and (f) function
fit performance.

The NNs used contain two hidden layers with 20 neurons, while other units are pre-
sented in table 2.



S. A.Shah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6481 17 of 19

Unite Initial value stopped value target value

Epoch 0 68 1000
Elapsed time - 00:00:00 -
Performance 28.9 5.26e− 12 0
Gradient 46.9 1.4e− 06 1e− 07
standard deviation 0.001 1e− 011 1e+ 10
Validation checks 0 0 6

Table 3: Training progress presentation for class V .

7. Conclusion

A discrete version of hepatitis B and C mathematical model was investigated by using
CFDO. Sufficient conditions were developed for the existence theory and stability analysis.
Moreover, numerical illustrations for various densities including healthy cells, infected and
virus free cells were presented against two different sets of fractional orders. Also some
stability results for the corresponding equilibrium points were given. In addition, NNs
scheme was established and different data dependence results by using the generated data
from the numerical scheme given in Eq. (20) using the provided initial values. Various
results related to MSE, RMSE, regression coefficient, and absolute errors for all three
classes were investigated. From the obtained analysis, it has been found that NNs can be
used as alternative powerful tools to investigate different epidemic diseases models. In the
future, we will extend our results to more complex dynamical systems of hepatitis B and
C mathematical models.

Conflict of interest

We have no conflict of interest.

Acknowledgements

The authors would like to acknowledge Prince Sultan University for paying the APC
and valuable support.

References

[1] Ching Shan Chou and Avner Friedman. Introduction to mathematical biology: mod-
eling, analysis, and simulations. Springer, 2016.

[2] MG Roberts and JAP Heesterbeek. Mathematical models in epidemiology, volume
215. EOLSS Abu Dhabi, United Arab Emirates, 2003.

[3] Shaher Momani, Ranbir Kumar, HM Srivastava, Sunil Kumar, and Samir Hadid. A
chaos study of fractional sir epidemic model of childhood diseases. Results in Physics,
27:104422, 2021.



S. A.Shah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6481 18 of 19

[4] A Isaacs and CH Andrewes. Spread of influenza. British medical journal, 2(4737):921,
1951.

[5] Ingrid A Binswanger, Komal J Narwaney, Jennifer C Barrow, Kathleen B Albers,
Laura Bechtel, Claudia A Steiner, Jo Ann Shoup, and Jason M Glanz. Association
between severe acute respiratory syndrome coronavirus 2 antibody status and rein-
fection: A case-control study nested in a colorado-based prospective cohort study.
Preventive Medicine Reports, 37:102530, 2024.

[6] A Melo-Florián. Spanish flu: learnings for xxith century. Medicine, pages 1–8, 2025.
[7] Helen S Te and Donald M Jensen. Epidemiology of hepatitis b and c viruses: a global

overview. Clinics in liver disease, 14(1):1–21, 2010.
[8] Don Ganem and Alfred M Prince. Hepatitis b virus infection -natural history and

clinical consequences. New England journal of medicine, 350(11):1118–1129, 2004.
[9] Syed Asad Ali, Rafe MJ Donahue, Huma Qureshi, and Sten H Vermund. Hepatitis

b and hepatitis c in pakistan: prevalence and risk factors. International journal of
infectious diseases, 13(1):9–19, 2009.

[10] Maud Lemoine and Mark Thursz. Hepatitis c, a global issue: access to care and new
therapeutic and preventive approaches in resource-constrained areas. In Seminars in
liver disease, volume 34, pages 089–097. Thieme Medical Publishers, 2014.

[11] Ngan Hoe Lee and Kit Ee Dawn Ng. Mathematical modelling: From theory to practice,
volume 8. World Scientific, 2015.

[12] Amit Huppert and Guy Katriel. Mathematical modelling and prediction in infectious
disease epidemiology. Clinical microbiology and infection, 19(11):999–1005, 2013.

[13] Rafael Pantoja Rangel, MDLGMagaña, Ricardo Ulloa Azpeitia, and Elena Nesterova.
Mathematical modeling in problem situations of daily life. Journal of Education and
Human Development, 5(1):62–76, 2016.

[14] Piotr Ostalczyk. Discrete fractional calculus: applications in control and image pro-
cessing, volume 4. World scientific, 2015.

[15] Abdulrahman Al-Khedhairi, Abdelalim A Elsadany, and Amr Elsonbaty. On the
dynamics of a discrete fractional-order cournot–bertrand competition duopoly game.
Mathematical Problems in Engineering, 2022(1):8249215, 2022.

[16] Amr Elsonbaty and AA Elsadany. On discrete fractional-order lotka-volterra model
based on the caputo difference discrete operator. Mathematical Sciences, 17(1):67–79,
2023.

[17] Lan-Lan Huang, Ju H Park, Guo-Cheng Wu, and Zhi-Wen Mo. Variable-order frac-
tional discrete-time recurrent neural networks. Journal of Computational and Applied
Mathematics, 370:112633, 2020.

[18] Zai-Yin He, Abderrahmane Abbes, Hadi Jahanshahi, Naif D Alotaibi, and Ye Wang.
Fractional-order discrete-time sir epidemic model with vaccination: Chaos and com-
plexity. Mathematics, 10(2):165, 2022.

[19] Vasily E Tarasov and Valentina V Tarasova. Long and short memory in economics:
fractional-order difference and differentiation. arXiv preprint arXiv:1612.07903, 2016.

[20] Yuexi Peng, Shaobo He, and Kehui Sun. Chaos in the discrete memristor-based
system with fractional-order difference. Results in Physics, 24:104106, 2021.



S. A.Shah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6481 19 of 19

[21] Fulai Chen and Yong Zhou. Existence and ulam stability of solutions for dis-
crete fractional boundary value problem. Discrete Dynamics in Nature and Society,
2013(1):459161, 2013.

[22] Pauline Van den Driessche. Reproduction numbers of infectious disease models. In-
fectious disease modelling, 2(3):288–303, 2017.

[23] Zahura Khatun, Md Shahidul Islam, and Uttam Ghosh. Mathematical modeling of
hepatitis b virus infection incorporating immune responses. Sensors International,
1:100017, 2020.

[24] Olga Igorevna Krivorotko, N Yu Zyatkov, and Sergey Igorevich Kabanikhin. Modeling
epidemics: Neural network based on data and sir-model. Computational Mathematics
and Mathematical Physics, 63(10):1929–1941, 2023.

[25] Behnam Nikparvar, Md Mokhlesur Rahman, Faizeh Hatami, and Jean-Claude Thill.
Spatio-temporal prediction of the covid-19 pandemic in us counties: modeling with a
deep lstm neural network. Scientific reports, 11(1):21715, 2021.

[26] Xiao Ning, Linlin Jia, Yongyue Wei, Xi-An Li, and Feng Chen. Epi-dnns: Epi-
demiological priors informed deep neural networks for modeling covid-19 dynamics.
Computers in biology and medicine, 158:106693, 2023.

[27] Muhammad Sher, Kamal Shah, Zeeshan Ali, Thabet Abdeljawad, and Manar
Alqudah. Using deep neural network in computational analysis of coupled systems
of fractional integro-differential equations. Journal of Computational and Applied
Mathematics, page 116912, 2025.

[28] Kamal Shah, Thabet Abdeljawad, Bahaaeldin Abdalla, and Zeeshan Ali. Analyzing a
coupled dynamical system of materials recycling in chemostat systems with artificial
deep neural network. Modeling Earth Systems and Environment, 11(5):1–14, 2025.

[29] Hasib Khan, Jehad Alzabut, and DK Almutairi. Applications of artificial intelligence
for clusters analysis of uranium decay via a fractional order discrete model. Partial
Differential Equations in Applied Mathematics, 13:101056, 2025.

[30] Hasib Khan, Jehad Alzabut, DK Almutairi, and Wafa Khalaf Alqurashi. The use of
artificial intelligence in data analysis with error recognitions in liver transplantation in
hiv-aids patients using modified abc fractional order operators. Fractal & Fractional,
9(1), 2025.