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

Criteria for Finite-Time Convergence in Discrete
Variable-Order Fractional FitzHugh–Nagumo

Reaction–Diffusion Systems
Nidal Anakira1,2,∗, Iqbal H. Jebril3, Iqbal M. Batiha3,4, Mohammad S. Hijazi5,∗,
Tala Sasa6

1 Mathematics Education Program, Faculty of Education and Arts, Sohar University, Sohar
311, Oman
2 Jadara University Research Center, Jadara University, Jordan
3 Department of Mathematics, Al Zaytoonah University of Jordan, Amman 11733, Jordan
4 Nonlinear Dynamics Research Center (NDRC), Ajman University, Ajman, UAE
5 Department of Mathematics, College of Sciences, Jouf University, Sakaka, Saudi Arabia
6 Department of Mathematics, Faculty of Science, Private Applied Science University,
Amman, Jordan

Abstract. This study addresses the problem of finite-time stability (FTS) for a discrete-time
FitzHugh–Nagumo reaction–diffusion system (FHN–RDs) governed by a variable-order (VO) Ca-
puto fractional difference operator. The discrete fractional formulation is obtained by combining
a central difference approximation for the spatial derivative with a VO fractional operator for the
temporal derivative. The analysis begins with proving the well-posedness of solutions for the pro-
posed discrete model. The main contribution lies in establishing an FTS criterion. By employing
a discrete fractional Gronwall-type inequality, we derive a sufficient stability condition expressed
through the discrete Mittag–Leffler function (MLF). Finally, a numerical simulation is provided
to illustrate the applicability of the theoretical findings, confirming that the system state remains
bounded within a prescribed limit over a finite time horizon.
2020 Mathematics Subject Classifications: 39A12, 34A08, 65M06, 92C15
Key Words and Phrases: FitzHugh–Nagumo system, finite-time stability, reaction-diffusion
systems, Gronwall’s inequality, Mittag-Leffler function

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

Email addresses: nanakira@su.edu.om (N. Anakira),
i.jebri@zuj.edu.jo (I. H. Jebril), i.batiha@zuj.edu.jo (I. M. Batiha), Mshijazi@ju.edu.sa (M. S.
Hijazi), t_sasa@asu.edu.jo (T. Sasa)

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



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7005 2 of 17

1. Introduction

Mathematical modeling plays a fundamental role in exploring and interpreting the
complex behavior of biological systems. Among the most prominent frameworks in this
domain are reaction–diffusion models, which capture how the spatial distribution of sub-
stances evolves under the combined influence of local reactions and diffusion. Turing’s
seminal work demonstrated that such systems can spontaneously produce patterns [1],
laying the foundation for morphogenesis and numerous spatio–temporal phenomena ob-
served in nature [2, 3]. A notable representative of this class is the FHN system, proposed
by FitzHugh [4] and independently developed by Nagumo and colleagues [5]. Conceived
as a reduction of the Hodgkin–Huxley model for nerve impulse transmission [6], the FHN
system retains essential features of excitable media, such as a stable equilibrium, an exci-
tation threshold, and a refractory phase [7, 8]. Owing to its simplicity and rich dynamics,
the FHN model has found applications far beyond neurophysiology, including wave propa-
gation and pattern formation in cardiology, chemical reactions, and population dynamics
[9, 10].

Conventional reaction–diffusion models typically rely on integer–order partial differ-
ential equations, inherently assuming that the underlying dynamics are Markovian, i.e.,
the future state depends solely on the current state. However, many real–world biological
and physical systems exhibit memory and hereditary characteristics, where the system’s
evolution also reflects its past history [11, 12]. Fractional calculus, which extends differ-
entiation and integration to arbitrary orders, offers a versatile framework for modeling
such non–local and memory–dependent processes [13–17]. Incorporating fractional deriva-
tives into diffusion models enables the description of anomalous diffusion [18–24], which
frequently arises in heterogeneous or crowded environments [25–27]. In reaction kinetics,
fractional operators can represent distributed delays and complex memory effects. Con-
sequently, FO-FHN models have been extensively investigated to study neuronal activity
with memory effects [28, 29].

To further increase flexibility, VO fractional calculus has been introduced [30, 31], in
which the derivative order varies with time, spatial position, or other parameters. This
allows for modeling systems whose memory characteristics or internal processes change
dynamically, such as in media with evolving properties [32, 33]. VO models are particularly
well–suited to capturing adaptive or evolving biological phenomena.

Since analytical solutions for such complex models are rarely attainable, numerical
techniques play a pivotal role in their analysis. Discrete fractional calculus provides a
rigorous foundation for approximating fractional operators and studying the corresponding
difference equations [34–36]. In this work, we formulate a discrete–time VO fractional
version of the FHN-RDs. A central aspect of system analysis is stability. For FO dynamics,
the classical notion of exponential stability is often replaced by Mittag–Leffler stability
[37–42]. In many engineering and biological contexts, however, the primary concern is not
the asymptotic state but ensuring that system trajectories remain within certain bounds
over a finite time horizon. This motivates the study of FTS [43], which has become an
important topic in fractional systems theory [44–46].



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7005 3 of 17

The present work focuses on establishing FTS conditions for the discrete VO fractional
FHN system. Using discrete fractional calculus tools, and in particular a discrete fractional
Gronwall–type inequality [47], we derive a sufficient criterion for stability. The remainder
of the paper is organized as follows: we first develop the discrete model and establish the
well-posedness of solutions; we then present the main FTS result; finally, computational
analysis are carried out to illustrate and validate the theoretical analysis. This paper is
organized as follows: Section 2 formulates the discrete model and establishes the well-
posedness of solutions. Section 3 presents the main FTS results. Section 4 provides
numerical simulations to illustrate and validate the theoretical analysis.

2. Problem Formulation

We consider the FHN–RDs, originally proposed in [48], in the following form:

∂u

∂t
= d1∆u− u3 + (β + 1)u2 − βu− v, x ∈ Ω, t > 0,

∂v

∂t
= d2∆v + ϵu− ϵγv, x ∈ Ω, t > 0,

∂xu = ∂xv = 0, x ∈ ∂Ω, t > 0,

u(x, 0) = u0(x), v(x, 0) = v0(x), x ∈ Ω.

(1)

Here, the variables, parameters, and operators appearing in system (1) are summarized in
Table 1:

Symbol / Parameter Description
Ω Bounded domain with smooth boundary ∂Ω
∆ Laplacian operator
u Membrane potential at (x, t) ∈ Ω× (0,∞)
v Combination of potassium activation and sodium inactivation
β Positive constant, 0 < β < 1

2
ϵ Positive constant, ϵ ≪ 1
γ Positive constant

To incorporate memory effects, we employ the VO Caputo FO derivative, resulting in
the following VO time–fractional FHN–RDs:{

C
0 D

δ(t)
t u− d1∆u = −u3 + (β + 1)u2 − βu− v,

C
0 D

δ(t)
t v − d2∆v = ϵu− ϵγv,

(2)

Here, the variables and parameters specific to the fractional formulation are summarized
in Table 2:



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7005 4 of 17

Symbol / Parameter Description
0 < δ(t) ≤ 1 Fractional VO
C
0 D

δ(t)
t Caputo fractional derivative

d1, d2 Diffusion coefficients

Let x ∈ [0, L], with spatial discretization xi+1 = xi + ∆x, i = 0, . . . ,m. Using the
central difference approximation, the second spatial derivative is given by

∂2y(x, t)

∂x2
≈ yi−1(t)− 2yi(t) + yi+1(t)

∆2
x

, y ∈ {u, v}. (3)

and using the second-order difference operator [49]:

∆2yi−1 = yi−1 − 2 yi + yi+1. (4)

we can write:
∂2y(x, t)

∂x2
≈ ∆2yi−1(t)

∆2
x

. (5)

Applying these approximations, we obtain the discrete VO fractional FHN–RD model:
C
ℏ ∆

δ(t)
t0

ui(t) =
d1
∆2

x
∆2Eθ(t)[ui−1](t)− Eθ(t)[u3i ](t) + (β + 1)Eθ(t)[u2i ](t)

−βEθ(t)[ui](t)− Eθ(t)[vi](t),

C
ℏ ∆

δ(t)
t0

vi(t) =
d2
∆2

x
∆2Eθ(t)[vi−1](t) + ϵEθ(t)[ui](t)− ϵγEθ(t)[vi](t).

(6)

where Eθ(t)[yi](t) := yi
(
t+ θ(t)

)
, with PBCs:

yj+m(t) = yj(t), j = 0, 1. (7)

and ICs:

yi(0) = ϕj(xj), j = 1, 2. (8)

Definition 1 ([49]). The VO Caputo fractional difference operator is defined as:
C
ℏ ∆

δ(t)
a κ(t) = ℏ∆

−(n−δ(t))
a ∆n

ℏκ(t), 0 < δ(t) ≤ 1, (9)

where

ℏ∆
−δ(t)
a κ(t) =

ℏ
Γ(δ(t))

s=a
ℏ∑

t
ℏ−δ(t)

(t− σ(sℏ))(δ(t)−1)
ℏ κ(sℏ), (10)

with σ(sℏ) = (s+ 1)ℏ and t ∈ (ℏN)a+δ(t)ℏ. The generalized falling factorial is:

t
(δ(t))
ℏ = ℏδ(t)

Γ( tℏ + 1)

Γ( tℏ − δ(t) + 1)
. (11)



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7005 5 of 17

Lemma 1 ([49]). The following properties hold:{
ℏ∆

−δ(t)
a+(1−δ(t))ℏ

C
ℏ ∆

δ(t)
a κ(t) = κ(t)− κ(a),

C
ℏ ∆

δ(t)κ = 0, 0 < δ(t) ≤ 1.
(12)

Lemma 2. The nonlinear fractional partial difference system (6) admits a unique solution
given by

uni = ϕ1,i +
ℏδn

Γ(δn)

n∑
j=1

w
(δn)
n−j

[
d1
∆2

x

∆2u j
i−1 −

(
u j
i

)3
+ (β + 1)

(
u j
i

)2 − β u j
i − v j

i

]
,

vni = ϕ2,i +
ℏδn

Γ(δn)

n∑
j=1

w
(δn)
n−j

[
d2
∆2

x

∆2v j
i−1 + ϵ u j

i − ϵγ v j
i

]
,

(13)

for 1 ≤ i ≤ m and n > 0, where

w
(δn)
n−j :=

Γ
(
n− j + δn

)
Γ(n− j + 1)

.

Proof 1. From system (6), we write:
C
ℏ ∆

δn
t0 ui(t) =

d1
∆2

x

∆2ui−1(t+ ℏδn) + F (ui(t+ ℏδn), vi(t+ ℏδn)) ,

C
ℏ ∆

δn
t0 vi(t) =

d2
∆2

x

∆2vi−1(t+ ℏδn) +G(ui(t+ ℏδn), vi(t+ ℏδn)) ,
(14)

where

F (ui, vi) = −u3i + (β + 1)u2i − βui − vi, (15)
G(ui, vi) = ϵui − ϵγvi. (16)

Using Lemma 1, we apply the fractional sum operator:
ℏ∆

−δn

t0+(1−δn)ℏ
C
ℏ ∆

δn
t0 ui(t) = ℏ∆

−δn

t0+(1−δn)ℏ

[ d1
∆2

x

∆2ui−1(t+ ℏδn) + F (·)
]
,

ℏ∆
−δn

t0+(1−δn)ℏ
C
ℏ ∆

δn
t0 vi(t) = ℏ∆

−δn

t0+(1−δn)ℏ

[ d2
∆2

x

∆2vi−1(t+ ℏδn) +G(·)
]
,

(17)

where (·) = (ui(t+ ℏδn), vi(t+ ℏδn)). By Definition 1 and Lemma 1, we obtain:
ui(t)− ϕ1,i =

ℏ
Γ(δn)

t
ℏ−δn∑

s=
t0
ℏ +1−δn

(t− σ(sℏ))δ
n−1

ℏ

[ d1
∆2

x

∆2ui−1(sℏ+ ℏδn) + F (·)
]
,

vi(t)− ϕ2,i =
ℏ

Γ(δn)

t
ℏ−δn∑

s=
t0
ℏ +1−δn

(t− σ(sℏ))δ
n−1

ℏ

[ d2
∆2

x

∆2vi−1(sℏ+ ℏδn) +G(·)
]
,

(18)



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7005 6 of 17

where (t− σ(sℏ))δ
n−1

ℏ = ℏδn−1
Γ
(
t
ℏ − s

)
Γ
(
t
ℏ − s− δn + 1

) . Substituting this yields:


ui(t) = ϕ1,i +

ℏδn

Γ(δn)

t
ℏ−δn∑

s=
t0
ℏ +1−δn

Γ
(
t
ℏ − s

)
Γ
(
t
ℏ − s− δn + 1

) [ d1
∆2

x

∆2ui−1(sℏ+ ℏδn) + F (·)
]
,

vi(t) = ϕ2,i +
ℏδn

Γ(δn)

t
ℏ−δn∑

s=
t0
ℏ +1−δn

Γ
(
t
ℏ − s

)
Γ
(
t
ℏ − s− δn + 1

) [ d2
∆2

x

∆2vi−1(sℏ+ ℏδn) +G(·)
]
,

(19)

where F (·) = F (ui(sℏ+ ℏδn), vi(sℏ+ ℏδn)) and similarly for G. Then,

ui(t) = ϕ1,i +
ℏδn

Γ(δn)

t
ℏ−δn∑

s=
t0
ℏ +1−δn

Γ
(
t
ℏ − s

)
Γ
(
t
ℏ − s− δn + 1

)
×
[ d1
∆2

x

(
ui−1(sℏ+ ℏδn)− 2ui(sℏ+ ℏδn) + ui+1(sℏ+ ℏδn)

)
−u3i (sℏ+ ℏδn) + (β + 1)u2i (sℏ+ ℏδn)− βui(sℏ+ ℏδn)− vi(sℏ+ ℏδn)

]
,

vi(t) = ϕ2,i +
ℏδn

Γ(δn)

t
ℏ−δn∑

s=
t0
ℏ +1−δn

Γ
(
t
ℏ − s

)
Γ
(
t
ℏ − s− δn + 1

)
×
[ d2
∆2

x

(
vi−1(sℏ+ ℏδn)− 2vi(sℏ+ ℏδn) + vi+1(sℏ+ ℏδn)

)
+ϵui(sℏ+ ℏδn)− ϵγvi(sℏ+ ℏδn)

]
.

(20)

The solution becomes:
uni = ϕ1,i +

ℏδn

Γ(δn)

n∑
j=1

w
(δn)
n−j

[ d1
∆2

x

∆2u j
i−1 −

(
u j
i

)3
+ (β + 1)

(
u j
i

)2 − β u j
i − v j

i

]
,

vni = ϕ2,i +
ℏδn

Γ(δn)

n∑
j=1

w
(δn)
n−j

[ d2
∆2

x

∆2v j
i−1 + ϵ u j

i − ϵγ v j
i

]
.

(21)

3. Main Results

Lemma 3. The linear fractional difference equation

C
ℏ ∆

δ(t)
t0

κ(t) = λκ
(
t+ δ(t)ℏ

)
,


0 < δ(t) ≤ 1,

|λ| < 1,

t ∈ (ℏN)t0 ,
κ(t0) = κ0.

(22)



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7005 7 of 17

admits a unique solution given by the discrete MLF:

κ(t) = κ0Eδ(t)(λ, t) = κ0

∞∑
j=0

λj

(
t
ℏ −

(
t0
ℏ + 1

)
+ jδ(t)

)(jδ(t))
ℏ

Γ(jδ(t) + 1)
. (23)

Proof 2. The solution is derived using Picard iteration. Define the recurrence:

κk+1(t) = κ0 + λ ℏ∆
−δ(t)
t0+(1−δ(t))ℏκk(s+ δ(t)ℏ), t ∈ (ℏN)t0 . (24)

The first iterates are:

κ1(t) = κ0 + λκ0

(
t
ℏ −

(
t0
ℏ + 1

)
+ δ(t)

)(δ(t))
ℏ

Γ(δ(t) + 1)
,

κ2(t) = κ0 + λκ0

(
t
ℏ −

(
t0
ℏ + 1

)
+ δ(t)

)(δ(t))
ℏ

Γ(δ(t) + 1)
+ λ2κ0

(
t
ℏ −

(
t0
ℏ + 1

)
+ 2δ(t)

)(2δ(t))
ℏ

Γ(2δ(t) + 1)
,

...

κn(t) = κ0

n∑
k=0

λk

(
t
ℏ −

(
t0
ℏ + 1

)
+ kδ(t)

)(kδ(t))
ℏ

Γ(kδ(t) + 1)
.

To prove convergence, consider the general term:

ak = λk

(
t
ℏ −

(
t0
ℏ + 1

)
+ kδ(t)

)(kδ(t))
ℏ

Γ(kδ(t) + 1)
. (25)

Using the Beta function property

B(x, y) =
Γ(x)Γ(y)

Γ(x+ y)
,

we simplify:

ak =
λk

kδ(t)B
(
t−t0
ℏ , kδ(t)

) . (26)

Applying the d’Alembert ratio test and Stirling’s approximation:

B

(
t− t0
ℏ

, kδ(t)

)
∼ Γ

(
t− t0
ℏ

)
(kδ(t))−(t−t0)/ℏ, (27)

yields:
lim
k→∞

∣∣∣∣ak+1

ak

∣∣∣∣ = |λ| < 1. (28)

Thus, the series converges absolutely to the solution (23).



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7005 8 of 17

Definition 2. System (6) is FTS with respect to {η, ε, J} (η < ε) if

∥ϕ1∥c + ∥ϕ2∥c < η (29)

implies
∥u(t)∥+ ∥v(t)∥ < ε, ∀t ∈ J = [t0, t0 + T ] ∩ (ℏN)t0 , (30)

where

∥ϕ∥c = max
1≤i≤m

|ϕ(xi)|, (31)

∥u(t)∥ =
m∑
i=1

|ui(t)|. (32)

Lemma 4. Let f(t), g(t) > 0, non-decreasing discrete functions, and assume g(t) ≤ M
for t ∈ J . If

κ(t) ≤ f(t) + g(t) ℏ∆
−δ(t)
t0+(1−δ(t))ℏκ(t+ δ(t)ℏ), t ∈ (ℏN)t0 , (33)

then
κ(t) ≤ f(t)Eδ(t)(g(t), t) , t ∈ (ℏN)t0 . (34)

Proof 3. Introduce the operator

Aφ(t) := g(t) ℏ∆
−δ(t)
t0+(1−δ(t))ℏφ(t+ δ(t)ℏ). (35)

Under this notation, (33) reads κ(t) ≤ f(t) +Aκ(t). By the monotonicity of A, it follows
that

κ(t) ≤
n−1∑
k=0

Akf(t) +Anκ(t). (36)

We claim that for every n ≥ 1,

Anκ(t) ≤ [g(t− (n− 1)δ(t)ℏ)]n ℏ∆
−nδ(t)
t0+(1−δ(t))ℏκ(t+ δ(t)ℏ). (37)

For n = 1, the inequality holds by definition. Assume it holds for some n = k. For
n = k + 1:

Ak+1κ(t) = A
(
Akκ(t)

)
≤ g(t− kδ(t)ℏ) ℏ∆

−δ(t)
t0+(1+(k−1)δ(t))ℏ

[
[g(t− (k − 1)δ(t)ℏ)]k ℏ∆

−kδ(t)
t0+(1−δ(t))ℏκ(t+ δ(t)ℏ)

]
= [g(t− kδ(t)ℏ)]k+1

ℏ∆
−(k+1)δ(t)
t0+(1−δ(t))ℏκ(t+ δ(t)ℏ),

by the composition property of fractional sums.
Since g(t) ≤ M and limn→∞Anκ(t) = 0, we deduce

Anf(t) ≤ f(t)
[g(t)]n

(
t
ℏ −

(
t0
ℏ + 1

)
+ nδ(t)

)(nδ(t))
ℏ

Γ(nδ(t) + 1)
,



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7005 9 of 17

κ(t) ≤
∞∑
k=0

Akf(t) ≤ f(t)Eδ(t)(g(t), t) .

Theorem 1. The system (6) is FTS if

Eδ(t)(ω, t) ≤
ε

η
, ∀t ∈ J, (38)

where ω is defined in (42).

Proof 4. From the solution representation and norm estimates:

∥u(t)∥ ≤ ∥ϕ1∥c + ℏ∆
−δ(t)
t0+(1−δ(t))ℏ

[(
β +

4d1
∆2

x

)
∥u(t+ δ(t)ℏ)∥+ ∥u(t+ δ(t)ℏ)∥3

+ (β + 1)∥u(t+ δ(t)ℏ)∥2 + ∥v(t+ δ(t)ℏ)∥
]
. (39)

Using the inequality ∥u∥3 + (β + 1)∥u∥2 ≤ ∥u∥ for ∥u∥ ≤ −(β+1)+
√

(β+1)2+4

2 :

∥u(t)∥ ≤ ∥ϕ1∥c + ℏ∆
−δ(t)
t0+(1−δ(t))ℏ

[(
4d1
∆2

x

+
β − 1 +

√
(β + 1)2 + 4

2

)
∥u(t+ δ(t)ℏ)∥

+ ∥v(t+ δ(t)ℏ)∥] . (40)

Similarly for v:

∥v(t)∥ ≤ ∥ϕ2∥c + ℏ∆
−δ(t)
t0+(1−δ(t))ℏ

[(
ϵγ +

4d2
∆2

x

)
∥v(t+ δ(t)ℏ)∥+ ϵ∥u(t+ δ(t)ℏ)∥

]
. (41)

Summing (40) and (41):

∥u(t)∥+ ∥v(t)∥ ≤ η + ℏ∆
−δ(t)
t0+(1−δ(t))ℏ

[(
ϵ+

4d1
∆2

x

+
β − 1 +

√
(β + 1)2 + 4

2

)
∥u(t+ δ(t)ℏ)∥

+

(
1 + ϵγ +

4d2
∆2

x

)
∥v(t+ δ(t)ℏ)∥

]
≤ η + ω ℏ∆

−δ(t)
t0+(1−δ(t))ℏ [∥u(t+ δ(t)ℏ)∥+ ∥v(t+ δ(t)ℏ)∥] ,

where

ω = max

{
ϵ+

4d1
∆2

x

+
β − 1 +

√
(β + 1)2 + 4

2
, 1 + ϵγ +

4d2
∆2

x

}
. (42)

Applying Lemma 4 yields:

∥u(t)∥+ ∥v(t)∥ ≤ ηEδ(t)(ω, t). (43)

Thus, Eδ(t)(ω, t) ≤ ε/η ensures FTS.



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7005 10 of 17

4. Numerical simulation

To validate the theoretical results established in Section 3, we present a comprehensive
numerical analysis of the FTS for the discrete VO fractional FHN system. The numer-
ical implementation follows the discrete framework developed in Section 2, with careful
attention to the accurate computation of the VO fractional operators and the spatial
discretization.

The discrete system in Eq. (6) was solved using an explicit numerical scheme based
on the Grünwald-Letnikov approximation for the VO fractional difference operator. The
spatial domain x ∈ [0, 10] was discretized with ∆x = 5 (yielding m = 2 spatial points due
to the periodic boundary conditions), and the temporal domain was discretized with time
step ∆t = 0.5 s. The VO fractional difference operator was computed using the definition
in Eq. (9), with the summation truncated at the current time step. For the numerical
evaluation of the discrete MLF in Eq.(??), we employed a truncated series representation
with 50 terms, which provided sufficient accuracy for the time intervals considered. The
norm calculations followed Definition 2, with ∥ · ∥c representing the maximum norm over
spatial points and ∥ · ∥ denoting the L1-norm as defined in Eq. (32).

We consider the following parameter values consistent with neuronal modeling appli-
cations:

(β, ε, γ, d1, d2,∆x) = (0.139, 0.45, 0.18, 0.5, 1, 5) (44)

With these parameters, the stability constant ω from Theorem 1 evaluates to:

ω = max

{
ε+

4d1
∆x2

+
β − 1 +

√
(β + 1)2 + 4

2
, 1 + εγ +

4d2
∆x2

}
= 1.25 (45)

The initial conditions were chosen as:{
ϕ1(xi) = 0.1(1 + sin(xi)),

ϕ2(xi) = 0.11(1 + sin(xi)),
(46)

which yield an initial norm:

∥ϕ1∥+ ∥ϕ2∥ = 0.096 < η (47)

where η = 0.1 represents the prescribed initial bound, and ε = 0.9 denotes the target
stability bound. Two distinct variable-order functions were examined to demonstrate the
flexibility of our stability framework:

• Case 1: δ1(t) = 0.3e−0.1t

• Case 2: δ2(t) =
0.2e−0.1t

t+ 1



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7005 11 of 17

For each case, we determined the maximum finite-time interval T for which the stability
condition in Eq. (38) holds:

Eδ1(t)(ω, t) ≤
ε

η
= 9 for t ∈ [0, T1] (48)

Eδ2(t)(ω, t) ≤
ε

η
= 9 for t ∈ [0, T2] (49)

Numerical evaluation of the discrete Mittag-Leffler function yielded:
T1 = 4.0 seconds with Eδ1(t)(ω, T1) = 8.355 (50)
T2 = 3.0 seconds withEδ2(t)(ω, T2) = 3.422 (51)

4.0.1. Results and Discussion

Figure 1 displays the evolution of ∥u(t)∥+∥v(t)∥ for Case 1 (δ1(t) = 0.3e−0.1t) over T1 = 4.0
seconds. The solution norm remains strictly below the stability bound ε = 0.9 throughout
the interval, confirming the finite-time stability as predicted by Theorem 1. The norm
initially increases to a peak value of approximately 0.85 at t = 2.5 seconds before gradually
decreasing, demonstrating the memory-dependent dynamics characteristic of fractional
systems.

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

Time, t (s)

0.0

0.2

0.4

0.6

0.8

1.0

N
or
m

V
al
ue

δ1(τ) = 0.3e−0.1τ

‖u(t)‖+ ‖v(t)‖ (Numerical)

ηEδ(τ)(ω, t) (Theoretical)

Stability Bound ε = 0.9

Figure 1: Estimation ∥u(t)∥+ ∥v(t)∥ within T = 4 s: δ1(t) = 0.3e−0.1t. The stability bound ε = 0.9 is shown
as a dashed red line.

Figure 2 presents the corresponding results for Case 2 (δ2(t) =
0.2e−0.1t

t+ 1
) over T2 =

3.0 seconds. Similar to Case 1, the solution norm stays within the prescribed bound



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7005 12 of 17

ε = 0.9 for the entire interval. However, the different variable-order function produces
a distinct dynamical profile, with a slightly higher peak value (approximately 0.88) at
t = 2.5 seconds. This difference highlights how the specific form of the variable-order
function influences the transient behavior while still preserving finite-time stability.

0.0 0.5 1.0 1.5 2.0 2.5 3.0

Time, t (s)

0.0

0.2

0.4

0.6

0.8

1.0

N
or
m

V
al
ue

δ2(τ) =
0.2e−0.1τ

τ + 1
‖u(t)‖+ ‖v(t)‖ (Numerical)

ηEδ(τ)(ω, t) (Theoretical)

Stability Bound ε = 0.9

Figure 2: Estimation ∥u(t)∥+ ∥v(t)∥ within T = 3 s: δ2(t) =
0.2e−0.1t

t+ 1
. The stability bound ε = 0.9 is shown

as a dashed red line.

The results validate our theoretical framework by demonstrating that:

(i) The system remains within the prescribed stability bound ε = 0.9 for the entire
finite-time interval [0, T ]

(ii) The actual finite-time interval T depends critically on the variable-order function
δ(t)

(iii) The discrete Mittag-Leffler function provides an accurate upper bound for the solu-
tion norm

(iv) Different variable-order functions produce different dynamical behaviors while main-
taining stability

To facilitate reproducibility of these results, Table 1 and Table 2 provide the numerical
values corresponding to Figures 1 and 2, respectively. These values were generated using
the Jupyter Notebook implementation described in Appendix A, which is available in the
supplementary materials.



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7005 13 of 17

Table 1: Numerical values for Figure 1: δ1(t) = 0.3e−0.1t, T = 4 seconds.

Time t (seconds) ∥u(t)∥+ ∥v(t)∥
0.0 0.2186
0.5 0.1952
1.0 0.1750
1.5 0.1604
2.0 0.1506
2.5 0.1454
3.0 0.1427
3.5 0.1414
4.0 0.1411

Table 2: Numerical values for Figure 2: δ2(t) =
0.2e−0.1t

t+ 1
, T = 3 seconds.

Time t (seconds) ∥u(t)∥+ ∥v(t)∥
0.0 0.2186
0.5 0.2080
1.0 0.2059
1.5 0.2064
2.0 0.2076
2.5 0.2088
3.0 0.2100

These findings underscore the practical utility of our stability criterion for predicting
the transient behavior of discrete fractional systems with variable-order dynamics. The
results also demonstrate that the VO framework offers enhanced modeling flexibility com-
pared to constant-order fractional models, as it can capture systems with evolving memory
characteristics. This is particularly relevant for biological applications where adaptation
and memory evolution are key features of the underlying processes.

5. Conclusion

In this work, we have presented a comprehensive investigation into the FTS of a dis-
crete FHN-RDs, distinguished by the incorporation of a VO Caputo fractional difference
operator. The study was motivated by the need for more sophisticated mathematical mod-
els that can accurately capture the complex memory effects and dynamic, time-varying
behaviors inherent in many biological and physical systems, such as neuronal networks.
The research commenced with the rigorous formulation of the discrete model. By em-
ploying a central difference scheme for the spatial Laplacian and a discrete VO fractional
operator for the temporal derivative, we systematically translated the continuous system
into a discrete-time framework suitable for numerical analysis. A crucial preliminary step



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7005 14 of 17

was to establish the well-posed of the solution for this newly formulated discrete system,
thereby ensuring a solid theoretical foundation for the subsequent stability analysis.

The principal contribution of this paper is the derivation of a novel, sufficient condition
for the FTS of the proposed model. This is particularly significant because, in many
practical applications, especially in biology and control engineering, guaranteeing that
system states remain within prescribed safe bounds over a finite operational interval is
of greater importance than assessing their long-term asymptotic behavior. By skillfully
applying a discrete fractional Gronwall-type inequality, a powerful tool in the analysis of
fractional difference equations, we established a stability criterion expressed elegantly in
terms of the discrete MLF. This result directly links the transient stability of the system
to its intrinsic memory properties, as captured by the fractional operator. Furthermore,
our analysis revealed a key insight: the interval of FTS is explicitly dependent on the
function defining the VO fractional order, δ(t). This highlights the profound impact
of dynamic memory on the system’s transient behavior and underscores the enhanced
modeling flexibility and descriptive power of the VO approach over traditional constant-
order fractional models. To validate our theoretical framework, a numerical example was
presented. This simulation not only corroborated the derived stability condition but also
provided a clear illustration of the system’s dynamics. By observing the norm of the
solution under different VO functions, we demonstrated that the system’s state can be
maintained within a predefined boundary over a finite time horizon, thereby confirming
the practical applicability of our results.

This study provides a rigorous and applicable framework for analyzing the transient
behavior of complex discrete systems with VO fractional dynamics. Looking ahead, this
work opens several avenues for future research. An immediate extension would be to inves-
tigate the influence of different classes of VO functions on system stability. Furthermore,
the introduction of stochastic perturbations could lead to an analysis of FTS in a noisy
environment, which is highly relevant for biological modeling. Finally, the conditions de-
rived herein could form the basis for designing control strategies aimed at enforcing FTS
in fractional-order systems. In summary, this research contributes valuable tools and in-
sights for the modeling and analysis of complex real-world phenomena where memory and
adaptation are key features.

References

[1] A. M. Turing. The chemical basis of morphogenesis. Philosophical Transactions of
the Royal Society of London. Series B, Biological Sciences, 237(641):37–72, 1952.

[2] J. D. Murray. Mathematical Biology II: Spatial Models and Biomedical Applications.
Springer, 2003.

[3] J. Keener and J. Sneyd. Mathematical Physiology. Springer, 1998.
[4] R. FitzHugh. Impulses and physiological states in theoretical models of nerve mem-

brane. Biophysical Journal, 1(6):445–466, 1961.
[5] J. Nagumo, S. Arimoto, and S. Yoshizawa. An active pulse transmission line simu-

lating nerve axon. Proceedings of the IRE, 50(10):2061–2070, 1962.



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7005 15 of 17

[6] A. L. Hodgkin and A. F. Huxley. A quantitative description of membrane current
and its application to conduction and excitation in nerve. The Journal of Physiology,
117(4):500–544, 1952.

[7] J. Rinzel. Excitation dynamics: insights from simplified membrane models. Federation
Proceedings, 40(13):2793–2804, 1981.

[8] C. Rocsoreanu, M. Sterpu, and A. Georgescu. The FitzHugh-Nagumo Model: Bifur-
cation and Dynamics. Springer, 2000.

[9] A. M. Pertsov, J. M. Davidenko, R. Salomonsz, W. T. Baxter, and J. Jalife. Spiral
waves of excitation underlie reentrant activity in isolated cardiac muscle. Circulation
Research, 72(3):631–650, 1993.

[10] D. E. Postnov, S. K. Han, T. G. Yim, and O. V. Sosnovtseva. Experimental obser-
vation of noise-induced coherence resonance in a an electrochemical system. Physical
Review E, 59(4):R3791, 1999.

[11] R. L. Magin. Fractional Calculus in Bioengineering. Begell House Publishers, 2006.
[12] R. Metzler and J. Klafter. The random walk’s guide to anomalous diffusion: a frac-

tional dynamics approach. Physics Reports, 339(1):1–77, 2000.
[13] I. Podlubny. Fractional Differential Equations. Academic Press, 1999.
[14] A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo. Theory and Applications of Frac-

tional Differential Equations. Elsevier, 2006.
[15] S. G. Samko, A. A. Kilbas, and O. I. Marichev. Fractional Integrals and Derivatives:

Theory and Applications. Gordon and Breach Science Publishers, 1993.
[16] I. Jebril, A. Lakehal, and S. Benyoussef. Fractional-order discrete predator–prey

system of leslie type: Existence, stability, and numerical simulation. International
Journal of Robotics and Control Systems, 5(2):1519–1538, 2025.

[17] Iqbal M Batiha, Hamzah O Al-Khawaldeh, Manal Almuzini, Waseem G Alshanti,
Nidal Anakira, and Ala Amourah. A fractional mathematical examination on breast
cancer progression for the healthcare system of jordan. Commun. Math. Biol. Neu-
rosci., 2025:Article–ID, 2025.

[18] I. M. Batiha, I. Bendib, A. Ouannas, P. Agarwal, N. Anakira, I. H. Jebril, and
S. Momani. Finite-time synchronization in a novel discrete fractional sir model for
covid-19. Special Issue on Advanced Computational Methods for Fractional Calculus,
43(1), 2025.

[19] A. Qazza, I. Bendib, R. Hatamleh, R. Saadeh, and A. Ouannas. Dynamics of the
gierer–meinhardt reaction–diffusion system: Insights into finite-time stability and
control strategies. Partial Differential Equations in Applied Mathematics, 14:Article
101142, 2025.

[20] Ahmed Bouchenak, Iqbal M Batiha, Mazin Aljazzazi, Nidal Anakira, Mohammad
Odeh, and Rasha Ibrahim Hajaj. Numerical and analytical investigations of frac-
tional self-adjoint equations and fractional sturm-liouville problems via modified con-
formable operator. Journal of Robotics and Control (JRC), 6(3):1410–1424, 2025.

[21] R. Hatamleh, I. Bendib, A. Qazza, R. Saadeh, A. Ouannas, and M. Dalah. Finite time
stability and synchronization of the glycolysis reaction-diffusion model. International
Journal of Neutrosophic Science, 25(4):371–386, 2025.



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7005 16 of 17

[22] S. Momani, I. M. Batiha, M. S. Hijazi, I. Bendib, A. Ouannas, and N. Anakira.
Fractional-order seir model for covid-19: Finite-time stability analysis and numerical
validation. International Journal of Neutrosophic Science, 26(1):266–282, 2025.

[23] H. Al-Taani, M. M. Hammad, O. Alomari, I. Bendib, and A. Ouannas. Finite-time
control of the discrete sel’kov–schnakenberg model: Synchronization and simulations.
AIP Advances, 15(2):025325, 2025.

[24] Khelifa Bouaziz, Nadhir Djeddi, Osama Ogilat, Iqbal M Batiha, Nidal Anakira, and
Tala Sasa. Stability analysis of a fractional-order lengyel–epstein chemical reaction
model. International Journal of Robotics and Control Systems, 5(2):1539–1551, 2025.

[25] I. M. Sokolov. Models of anomalous diffusion in crowded environments. Soft Matter,
8(35):9043–9052, 2012.

[26] A. N. Anber and Z. Dahmani. The laplace decomposition method for solving nonlin-
ear conformable fractional evolution equations. Int. J. Open Probl. Comput. Math.,
17(1):67–81, 2024.

[27] A. Anber and Z. Dahmani. The ldm and the cvim methods for solving time and
space fractional wu-zhang differential system. Int. J. Open Probl. Comput. Math.,
17(3):1–18, 2024.

[28] W. Teka, G. F. T. del Castillo-Negrete, and B. A. van Gorder. Anomalous transport
in the fractional fitzhugh-nagumo model. Communications in Nonlinear Science and
Numerical Simulation, 19(10):3539–3552, 2014.

[29] N. H. Sweilam, M. M. Khader, and R. F. Al-Bar. Numerical studies for a multi-term
fractional-order stencil for the fitzhugh-nagumo model. Journal of Advanced Research,
3(2):121–127, 2012.

[30] C. F. Lorenzo and T. T. Hartley. Variable order and distributed order fractional
operators. Nonlinear Dynamics, 29(1-4):57–98, 2002.

[31] C. F. M. Coimbra. Mechanics with variable-order differential operators. Annalen der
Physik, 12(11-12):692–703, 2003.

[32] H. Sun, Y. Zhang, D. Baleanu, W. Chen, and Y. Chen. A new collection of real
world applications of fractional calculus in science and engineering. Communications
in Nonlinear Science and Numerical Simulation, 64:213–231, 2018.

[33] L. E. S. Ramirez and C. F. M. Coimbra. On the variable order dynamics of the
nonlinear viscous-plastic friction in slender asperities. Journal of Sound and Vibration,
329(16):3277–3288, 2010.

[34] F. M. Atici and P. W. Eloe. A transform method in discrete fractional calculus.
International Journal of Difference Equations, 2(2):165–176, 2007.

[35] C. Goodrich and A. C. Peterson. Discrete Fractional Calculus. Springer, 2015.
[36] J. Cheng. Discrete Fractional Calculus. Springer, 2011.
[37] Y. Li, Y. Chen, and I. Podlubny. Mittag-leffler stability of fractional order nonlinear

dynamic systems. Automatica, 45(8):1965–1969, 2009.
[38] I. Stamova. Stability Analysis of Impulsive Functional Differential Equations. De

Gruyter, 2016.
[39] A. S. Hussain, K. D. Pati, A. K. Atiyah, and M. A. Tashtoush. Rate of occurrence

estimation in geometric processes with maxwell distribution: A comparative study



N. Anakira et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7005 17 of 17

between artificial intelligence and classical methods. Int. J. Adv. Soft Comput. Appl.,
17(1):1–15, 2025.

[40] M. Berir. Analysis of the effect of white noise on the halvorsen system of variable-
order fractional derivatives using a novel numerical method. Int. J. Adv. Soft Comput.
Appl., 16(3):294–306, 2024.

[41] P. Singh, N. Zade, P. Priyadarshi, and A. Gupte. The application of machine learning
and deep learning techniques for global energy utilization projection for ecologically
responsible energy management. Int. J. Adv. Soft Comput. Appl., 17(1):49–66, 2025.

[42] E. A. Mohammed and A. Lakizadeh. Benchmarking vision transformers for satellite
image classification based on data augmentation techniques. Int. J. Adv. Soft Comput.
Appl., 17(1):98–114, 2025.

[43] P. Dorato. Short-time stability in linear time-varying systems. In Proceedings of the
IRE International Convention Record, pages 83–87, 1961. Part 4.

[44] J. Shen, J. Cao, and S. M. Hayat. Finite-time stability of fractional-order complex-
valued neural networks with time delays. Neural Networks, 56:31–38, 2014.

[45] M. P. Lazarevic. Finite time stability analysis of fractional order time delay systems.
Masinstvo, 26(1):11–20, 2007.

[46] C. Fei, H. Li, and J. Yan. Finite-time stability of fractional-order neural networks
with delay. Neurocomputing, 182:161–166, 2016.

[47] T. Abdeljawad and D. Baleanu. On fractional derivatives with exponential kernel and
their discrete versions. Reports on Mathematical Physics, 84(3):367–378, 2019.

[48] I. M. Batiha, I. Bendib, A. Ouannas, I. H. Jebril, S. Alkhazaleh, and S. Momani. On
new results of stability and synchronization in finite-time for fitzhugh-nagumo model
using gronwall inequality and lyapunov function. Journal of Robotics and Control
(JRC), 5(6):1897–1909, Oct. 2024.

[49] I. Bendib, I. M. Batiha, A. Hioual, N. Anakira, M. Dalah, and A. Ouannas. On a
new version of gierer-meinhardt model using fractional discrete calculus. Results in
Nonlinear Analysis, 7(2):1–5, 2024.