Mathematics in Applied Sciences and Engineering https://ojs.lib.uwo.

ca/mase
Volume 4, Number 1, March 2023, pp.40-60 https://doi.org/10.5206/mase/15355

CHAOTIC DYNAMICS OF THE FRACTIONAL ORDER SCHNAKENBERG

MODEL AND ITS CONTROL

MD. JASIM UDDIN AND S. M. SOHEL RANA

Abstract. The Schnakenberg model is thought to be the Caputo fractional derivative. A discretiza-
tion process is first used to create caputo fractional differential equations for the Schnakenberg model.
The fixed points in the model are categorized topologically. Then, we show analytically that a frac-
tional order Schnakenberg model supports a Neimark-Sacker (NS) bifurcation and a Flip-bifurcation
under certain parametric conditions. Using central manifold and bifurcation theory, we demonstrate
the presence and direction of NS and Flip bifurcations. The parameter values and the initial conditions
have been found to profoundly impact the dynamical behavior of the fractional order Schnakenberg
model. Numerical simulations demonstrate chaotic behaviors like bifurcations, phase portraits, period
2, 4, 7, 8, 10, 16, 20 and 40 orbits, invariant closed cycles, and attractive chaotic sets in addition to val-
idating analytical conclusions. To support the system’s chaotic characteristics, we also quantitatively
compute the maximal Lyapunov exponents and fractal dimensions. Finally, the chaotic trajectory of
the system is stopped using the OGY approach, hybrid control method, and state feedback method.

1. Introduction
Differentiation and integration to arbitrary order, commonly known as fractional calculus, has re-
ceived a lot of interest from researchers. A mathematical notion from the 17th century is fractional calcu-
lus. But it may be regarded as a new research subject. Due to their close resemblance to memory-based
systems, which are present in most biological systems, fractional-order differential equations (FD) are
the most often employed [13]. It is possible to successfully explain fractional-order differential equations
in several fields, including science, engineering, finance, economics, and epidemiology[16, 17, 18, 19, 30].
Switching from an integer-order model to a fractional-order one necessitates accuracy in the order of
differentiation; even a slight change in α might greatly influence the final result[7]. Fractional Differen-
tial Equations can describe phenomena that IDEs can’t wholly model [20]. The complicated dynamics
of chaos and bifurcation may be seen in a nonlinear fractional differential system, much like in a non-
linear differential system. It’s fascinating and engaging to study disorder in fractional-order dynamical
systems[1, 4, 8, 9, 13].There are several strategies to apply the concept of differentiation to arbitrary or-
der. The most popular definitions are those by Riemann-Liouville, Caputo, and Grünwald-Letnikov[35].
Along with these definitions, researchers constantly look for the most effective strategy when developing
or altering their models, including specific numerical approaches[6, 21, 31].
Numerous discrete systems have aroused the interest of academics investigating the Neimark-Sacker
and flip bifurcations, stable orbits, and chaotic attractors (see [24, 25, 36, 37]). The center manifold
theory and standard form can mathematically quantify these phenomena.

Received by the editors 3 October 2022; accepted 3 March 2023; published online 22 March 2023.
2020 Mathematics Subject Classification. 37C25, 37D45, 39A28, 39A33.
Key words and phrases. Fractional order Schnakenberg model, flip and Neimark-Sacker bifurcations, maximum Lya-
punov exponent, fractal dimension, chaos control.

40
 CHAOTIC DYNAMICS OF THE FRACTIONAL ORDER SCHNAKENBERG MODEL AND ITS CONTROL 41

There are numerous cyclical and oscillatory processes in nature. Circadian rhythms, present in almost
every aspect of life, are possibly the most well-known of these recurring occurrences. A device made in
1979 by J. Schnakenberg demonstrated sustained oscillations for a straightforward glycolysis model (a
metabolic process that transforms glucose into cellular energy) [38], remarkably similar to a system of
four reactions called the ”Brusselator”[11]. A semi-analytical approach (see [5]) was used to investigate
the Schnakenberg model of a reaction-diffusion cell.
A chemical Schnakenberg model was investigated [23], which depicts autochemical processes with
rhythmic behavior that may have various biological and biochemical applications. A variable-order
space-time fractional reaction-diffusion Schnakenberg model’s numerical solutions were studied in [15].
In addition,the Schnakenberg model was described in [23, 29, 28], where a variety of numerical methods
were employed to get approximations of solutions.
The Schnakenberg model is an oscillatory chemical reaction model that Schnakenberg developed from
a number of hypothetical tri-molecular autocatalytic processes. In order to find the fewest possible
reactions and reactants that display limit-cycle behavior, this model has been developed. Schnakenberg
established that this kind of model needs at least three reactions, including one auto-catalytic reaction.
Thus, the following reaction scheme is developed for generic chemicals A, X, Y , and B:

m1 m m
X A, B →2 Y, 2X + Y →3 3X,
m−1

where X, Y represent chemicals with different concentrations and A, B have constant concentrations.
The following differential equations are implied by mass action theory

dNX 2
= m1 NA − m−1 NX + m3 NX NY ,
dτ (1.1)
dNY 2
= m2 NB − m3 NX NY ,
dτ
where NA , NB , NX and NY represent numbers of molecules for fixed and varying concentrations A, B, X
and Y respectively.
The non-dimensional form of hypothetical Schnackenberg Model comes from the above system.

ẋ = a − x + x2 y,
(1.2)
ẏ = b − x2 y,

where
r r
m3 m3
x= NX , y = NY , τ −→ t = τ m−1 ,
m−1 m−1
r r
m1 m3 m2
a= NA , b = NB .
m−1 m−1 m−1

Due to their effective computational outcomes and complex dynamical behavior, discrete-time models
governed by difference equations are preferable to continuous ones[32, 2]. Additionally, this reasoning
holds true for nonlinear oscillatory behavior associated with chemical reactions[22, 34, 14]. As a result,
we research the stability, chaos management, and bifurcation analysis of discrete counterparts of(1.2) .
Recently, a few authors[1, 4, 8, 9, 13, 3, 12] use the well-known Caputo fractional derivative instead of
ordinary derivatives in continuous model. Instead of using regular derivatives, the Schnakenberg model
can employ fractional derivatives because it works the same way. In one sense, this may interpret that
the rate of change for the concentrations of the chemical products will be slower and this may lead a
42 MD. JASIM UDDIN AND S. M. SOHEL RANA

better mathematical approximation. The fractional order Schnackenberg model is given as follows
Dα x(t) = a − x(t) + x2 (t)y(t),
(1.3)
Dα y(t) = b − x2 (t)y(t),
where α is the fractional order that satisfies α ∈ (0, 1] and t > 0. There are a lot of approaches for
discretize such kind of system. The piecewise constant approximation[3, 12] is one of them. The model
is discretized by using this method. The steps are listed below:
Let system (1.2) ’s starting conditions be x(0) = x0 , y(0) = y0 . The discretized version of system
(1.3) is given as:

t t t
Dα x(t) = a − x([ ]ρ) + x2 ([ ]ρ)y([ ]ρ),
ρ ρ ρ
(1.4)
α 2 t t
D y(t) = b − x ([ ]ρ)y([ ]ρ).
ρ ρ
t
First, let t ∈ [0, ρ), so ρ ∈ [0, 1). Thus, we obtain

Dα x(t) = a − x0 + x20 y0 ;
(1.5)
Dα y(t) = b − x20 y0 .
The solution of (1.5) is simplified to

tα
x1 (t) = x0 + J α a − x0 + x20 y0 = x0 + a − x0 + x20 y0 ,



αΓ(α)
α (1.6)
α t
b − x20 y0 = y0 + b − x20 y0 .



y1 (t) = y0 + J
αΓ(α)
t
Second, let t ∈ [ρ, 2ρ), so ρ ∈ [1, 2). Then

Dα x(t) = a − x1 + x21 y1
(1.7)
Dα y(t) = b − x21 y1
which have the following solution

x2 (t) = x1 (ρ) + Jρα a − x1 + x21 y1

(t − ρ)α
a − x1 + x21 y1 ,


= x1 (ρ) +
αΓ(α)
(1.8)
y2 (t) = y1 (ρ) + Jρα b − x21 y1

(t − ρ)α
b − x21 y1 ,


= y1 (ρ) +
αΓ(α)
where Z t
1
Jρα ≡ (t − τ )α−1 dτ, α > 0.
Γ(α) ρ
The result of doing the discretization process n times over

(t − nρ)α
a − xn (nρ) + x2n (nρ)yn (nρ) ,


xn+1 (t) = xn (nρ) +
αΓ(α)
(1.9)
(t − nρ)α
b − x2n (nρ)yn (nρ) ,


yn+1 (t) = yn (nρ) +
αΓ(α)
where t ∈ [nρ, (n + 1)ρ). For t −→ (n + 1)ρ, system (1.9) is reduced to
 CHAOTIC DYNAMICS OF THE FRACTIONAL ORDER SCHNAKENBERG MODEL AND ITS CONTROL 43

ρα
a − xn + x2n yn ,


xn+1 = xn +
Γ(α + 1)
(1.10)
ρα
b − x2n yn .


yn+1 = yn +
Γ(α + 1)
An ecological system in the real world is not always stable. A small change of control parameter
may destabilize the system from locally stable coexistence to producing chaotic orbits. In a discrete
dynamical system, Flip bifurcation and Neimark-Sacker bifurcation are the important mechanisms for
the generation of complex dynamics. Because in the discrete predator-prey system, both bifurcations
cause the system to jump from stable window to chaotic states through periodic and quasi-periodic
states, and trigger a route to chaos. When a system undergoes a Flip bifurcation, a sequence of
period-doubling cascades leads the system from steady state to chaos. On the other hand, when system
undergoes Neimark-Sacker bifurcation, it instigates a route to chaos, through a dynamic transition from
a stable state, to invariant closed cycle, with periodic and quasi-periodic states occurring in between,
to chaotic sets.
This paper’s remaining sections are organized as follows.: Sect. 2 investigates the fixed point topo-
logical classifications. In Sect. 3, we explore analytically the possibility that the system (1.10) would
experience a Flip or NS bifurcation under a certain parametric condition. In Sect.4, we numerically
show system dynamics that includes bifurcation diagrams, phase portraits, and MLEs to support our
analytical conclusions.To stabilize the chaos of the unmanaged system, we employ the OGY approach,
hybrid control method, and state feedback technique in Sect. 5. Sect. 6 presents a brief discussion.

2. Stability of Fixed point

The system (1.10) has a unique fixed point E(x∗ , y ∗ ), where x∗ = a + b and y ∗ = b/(a + b)2 , which
always exists for all permissible parameter values.
System (1.10)’s Jacobian matrix, evaluated at E(x∗ , y ∗ ), is as follows:
  α
 
ρ 2 ρα
1 + (−1 + 2x∗ y ∗ ) Γ(α+1) −x∗ Γ(α+1)
J (x∗ , y ∗ ) =  ρα
 2
 . (2.1)
−2x∗ y ∗ Γ(α+1) 1 − x∗
Now at E(x∗ , y ∗ ), the Jacobian matrix is given by

(−a+b) ρα ρ α
(a + b)2 Γ(α+1)
!
1+ a+b Γ(α+1)
JE = 2b ρα

ρα
 . (2.2)
− a+b Γ(α+1) 1 − (a + b)2 Γ(α+1)
The characteristic polynomial of the Jacobian Matrix can be written as

F̃ (λ) := λ2 − Tr(JE ) λ + Det(JE ) = 0, (2.3)

where T r(JE ) is the trace and Det(JE ) is the determinant of the Jacobian matrix JE , and is given by

(−a + b) ρα
Tr(JE ) =2 + − (a + b)2 ,
a+b Γ(α + 1)
(2.4)
a3 (−1 + ρb)b
ρ + 3a2 b(−1 + ρb)b
ρ + a(−1 + ρb)(−1 + 3b2 ρb) + b(1 + ρb − b2 ρb + b2 ρb2 )
Det(JE ) = .
a+b
where ρb = ρα /Γ(α + 1).
44 MD. JASIM UDDIN AND S. M. SOHEL RANA

The eigenvalues of the system(2.3) can be written as

p
T r(JE ) ± (T r(JE ))2 − 4Det(JE )
λ1,2 = .
2
By the Jury Criterion, the stability condition for the equlibrium point E(x∗ , y ∗ ) is given as follows:
F̃ (1) > 0, F̃ (−1) > 0, F̃ (0) − 1 < 0.
Let
√ ! α1
 
 −A2e ± L 
P DE = (a, b, ρ, α) : ρ = Γ(1 + α) = ρ± , L ≥ 0
 A1e 
where
A1e = (a + b)2 ,
a + a3 − b + 3a2 b + 3ab2 + b3 a − b + (a + b)3 A4e
A2e = − =− =− ,
a+b a+b a+b
A3e = 4
A4e = a − b + (a + b)3 ,
L = A22e − A1e A3e .
The system (1.10), however, undergoes a flip bifurcation at E when the parameters (a, b, ρ, α) fluctuate
within a narrow region of P DE .
Also let (  1 )
−A2e α
N SE = (a, b, ρ, α) : ρ = Γ(1 + α) = ρN S , L < 0 .
A1e
If the parameters (a, b, ρ, α) vary around the set N SE , system (1.10) will suffer an NS bifurcation at
that point.

Lemma 2.1. For every parameter value selection, the fixed point E is a
− sink if (i)L ≥ 0, ρ < ρ− (stable node),
(ii)L < 0, ρ < ρN S (stable focus),
− source (i)L ≥ 0, ρ > ρ+ (unstable node),
(ii)L < 0, ρ > ρN S (unstable focus),
− non-hyperbolic (i)L ≥ 0, ρ = ρ− or ρ = ρ+ ( saddle with flip),
(ii)(i)L < 0, ρ = ρN S ( focus),
− saddle: otherwise

3. Bifurcation Analysis
In this section, we will study the presence, direction, and stability analysis of flip and NS bifurcations
near to the fixed point E using center-manifold and bifurcation theory[26, 39, 40]. In other words, we
take ρ to be the bifurcation parameter.

3.1. Flip Bifurcation. We arbitrarily select the parameters (a, b, ρ, α) to locate in P DE . Consider
the system (1.10) at equilibrium point E(x∗ , y ∗ ). Assume the parameters lie in P DE .
Let L ≥ 0 and
√ ! α1
−A2e − L
ρ = ρ− = Γ(1 + α) .
A1e
Then, the eigenvalues of J(E) are given as
λ1 = −1 and λ2 = 3 + A2 ρ− .
 CHAOTIC DYNAMICS OF THE FRACTIONAL ORDER SCHNAKENBERG MODEL AND ITS CONTROL 45

Also note that the condition |λ2 6= 1| is equivalent to

A2 ρ− 6= −2, −4. (3.1)

Next, we use the transformation x̂ = x − x+ , ŷ = y − y + and set A(ρ) = J(x∗ , y ∗ ). We shift the fixed
point of system(1.10) to the origin. So the system(1.10) can be written as

     
x̂ x̂ F1 (x̂, ŷ, ρ− )
→ A(ρ− ) + (3.2)
ŷ ŷ F2 (x̂, ŷ, ρ− )

where X = (x̂, ŷ)T and

1 h √ i
a − b + (a + b)3 − L x̂ 2b3 ŷ + a2 (2a + x̂)ŷ + b2 (6a + x̂)ŷ


F1 (x̂, ŷ, ρ− ) = 5
(a + b)
1 h
3
√ i
L x̂ b(x̂ + 6a2 ŷ + 2ax̂ŷ) ,


+ 5
a − b + (a + b) −
(a + b)
√ i (3.3)
1 h
3
L x̂ 2b3 ŷ + a2 (2a + x̂)ŷ + b2 (6a + x̂)ŷ


F2 (x̂, ŷ, ρ− ) = − a − b + (a + b) −
(a + b)5
1 h √ i
a − b + (a + b)3 − L x̂ b(x̂ + 6a2 ŷ + 2ax̂ŷ) .


− 5
(a + b)

The system(3.2) can be expressed as

1 1 
4

Xn+1 = AXn + B (Xn , Xn ) + C (Xn , Xn , Xn ) + O kXn k
2 6

where
   
B1 (x, y) C1 (x, y, u)
B(x, y) = and C(x, y, u) =
B2 (x, y) C2 (x, y, u)

are symmetric multi-linear vector functions of x, y, u ∈ R2 and defined as follows:

2
X δ 2 F1 (ξ, ρ) 1 h
3
√ i 3

B1 (x, y) = xj yk = a − b + (a + b) − L a (x2 y1 + x1 y2 )
δξj δξk (a + b)5
j,k=1
ξ=0
1
3a2 b(x2 y1 + x1 y2 )


+ 5
(a + b)
1
b3 (x2 y1 + x1 y2 ) + 3a2 b(x2 y1 + x1 y2 ) + bx1 y1 ,


+ 5
(a + b)
2
X δ 2 F2 (ξ, ρ) 1 h √ i
a − b + (a + b)3 − L (a + b)3 (x2 y1 + x1 y2 ) + bx1 y1 ,


B2 (x, y) = xj yk = 5
δξj δξk (a + b)
j,k=1
ξ=0
46 MD. JASIM UDDIN AND S. M. SOHEL RANA

and
2
X δ 2 F1 (ξ, ρ) 1  √ 
C1 (x, y, u) = xj yk ul = 3
a − b + (a + b)3 − L (u1 x1 y2 )
δξj δξk δξl (a + b)
j,k,l=1
ξ=0
1 
3
√ 
+ a − b + (a + b) − L (u1 x1 y2 ) ,
(a + b)3
2
X δ 2 F1 (ξ, ρ) 1  √ 
C2 (x, y, u) = xj yk ul = − 3
a − b + (a + b)3 − L (u2 x1 y1 + u1 x2 y1 )
δξj δξk δξl (a + b)
j,k,l=1
ξ=0
1 
3
√ 
− a − b + (a + b) − L (u1 x1 y2 ) .
(a + b)3
Let q1 , q2 ∈ R2 be two eigenvectors of A and AT for eigenvalue λ1 (ρ− ) = −1 such that A (ρ− ) q1 = −q1
and AT (ρ− ) q2 = −q2 . Then by direct calculation we get
3 3
√ !
(−a−3b−+(a+b) −
− (a+b)2b(a−b+(a+b) L)  
√
3 − L)
q11
q1 = = ,
1 1
3
√ !
− a+3b−(a+b) + L  
√
a−b+(a+b)3 − L
q21
q2 = = .
1 1
In order to get hq1 , q2 i = 1, where hq1 , q2 i = q11 q21 + q12 q22 , we have to use the normalized vector
q2 = γF q2 , with γF = 1/(1 + q11 p11 ).
To obtain the direction of the flip bifurcation, we have to check the sign of l1 (ρ− ), the coefficient of
critical normal form ([26]) and is given by
1 1
q2 , B q2 , (A − I)−1 B(q1 , q1 ) .


l1 (ρ− ) = hq2 , C(q1 , q1 , q1 )i − (3.4)
6 2
In light of the reasoning above, the following theorem may be used to demonstrate the direction and
stability of Flip bifurcation.

Theorem 3.1. Suppose (3.1) holds and l1 (ρ− ) 6= 0, then Flip bifurcation at fixed point E(x∗ , y ∗ )
for system (2.1) if the ρ changes its value in small neighbourhood of P DE . Moreover, if l1 (ρ− ) <
0 (resp. l1 (ρ− ) > 0), then there is a smooth closed invariant curve that can bifurcate from E(x∗ , y ∗ ),
and the bifurcation is sub-critical (resp. super-critical).

3.2. Neimark-Sacker Bifurcation. When the parameters (a, b, α, ρ) ∈ N SE and L < 0, then the
eigenvalues of system (2.3) are given by
p
T r(JE ) ± i 4Det(JE ) − T r(JE )2
λ, λ̄ = (3.5)
2
where
(a − b)(a − b + (a + b)3 ) a − b + (a + b)3
T r(JE ) = 2 − − , Det(JE ) = 1.
(a + b)4 a+b
Let ρN S = −A2e /A1e . Then, the transversality and non-resonance conditions respectively read
d|λi (ρ)|
dρ |ρ=ρN S = − A22e 6= 0;
A22 (3.6)
−(T r(JE ))|ρ=ρN S 6= 0 ⇒ A1e 6= 2, 3.
 CHAOTIC DYNAMICS OF THE FRACTIONAL ORDER SCHNAKENBERG MODEL AND ITS CONTROL 47

Using the transformation x̂ = x − x+ , ŷ = y − y + and setting A(ρ) = J(x∗ , y ∗ ), we move the system(
1.10)’s fixed point to the origin. So the system(1.10) can be written as

X → A(ρ)X + F (3.7)
where X = (x̂, ŷ)T and F = (F1 , F2 )T are given by

(a − b + (a + b)3 )x̂ 2b3 ŷ + a2 (2a + x̂)ŷ + b2 (6a + x̂)ŷ + b(x̂ + 6a2 ŷ + 2ax̂ŷ)
F1 (x̂, ŷ, ρN S ) = ,
(a + b)5

(3.8)
(a − b + (a + b)3 )x̂ 2b3 ŷ + a2 (2a + x̂)ŷ + b2 (6a + x̂)ŷ + b(x̂ + 6a2 ŷ + 2ax̂ŷ)
F2 (x̂, ŷ, ρN S ) = − .
(a + b)5
The system(3.2) can be expressed as

1 1 
4

Xn+1 = AXn + B (Xn , Xn ) + C (Xn , Xn , Xn ) + O kXn k
2 6
where
2(a − b + (a + b)3 ) 
(a + b)3 (x2 y1 + x1 y2 ) + bx1 y1 ,

B1 (x, y) = 5
(a + b)
2(a − b + (a + b)3 ) 
(a + b)3 (x2 y1 + x1 y2 ) + bx1 y1 ,

B2 (x, y) = − 5
(a + b)
and
2(a − b + (a + b)3 )
C1 (x, y, u) = (u2 x1 y1 + u1 x2 y1 + u1 x1 y2 ) ,
(a + b)5
2(a − b + (a + b)3 )
C2 (x, y, u) = − (u2 x1 y1 + u1 x2 y1 + u1 x1 y2 ) .
(a + b)5
Let q1 , q2 ∈ C2 be two eigenvectors of A(ρN S ) and AT (ρN S )for eigenvalue λ(ρN S ) and λ̄(ρN S ) such
that

A (ρNS ) q1 = λ (ρNS ) q1 , A (ρNS ) q¯1 = λ̄ (ρNS ) q¯1 ,

(3.9)
AT (ρNS ) q2 = λ̄ (ρNS ) q2 , AT (ρNS ) q¯2 = λ (ρNS ) q¯2 .
Therefore, using a straight calculation, we obtain
3
√ !
− −a+b+(a+b) + L
 
ζ1 + iζ2
q1 = 4b = ,
1 1
√ !
−a+b+(a+b)3 − L  
2(a+b)3
ξ1 + iξ2
q2 = = .
1 1
where √
−a + b + (a + b)3 − −L
ζ1 = − , ζ2 = ,
4b 4b
√
−a + b + (a + b)3 − −L
ξ1 = 3
, ξ2 = .
2(a + b) 2(a + b)3
For hq1 , q2 i = 1 to hold, we can take normalized vector as q2 = γN S q2 where
1
γN S = .
1 + (ξ1 + iξ2 )(ζ1 − iζ2 )
The following is how the eigenvectors are calculated:
48 MD. JASIM UDDIN AND S. M. SOHEL RANA

 
ζ1 + iζ2
q1 = ,
1
ξ1 +iξ2
!
1+(ξ1 +iξ2 )(ζ1 −iζ2 )
q2 = 1 .
1+(ξ1 +iξ2 )(ζ1 −iζ2 )

We break down X ∈ R2 as X = zq1 + z̄ q¯1 by considering ρ vary near to ρN S and for z ∈ C. The
precise formulation of z is z = hq2 , Xi. As a result, the system (2.3) changed to the following system
for |ρ| near ρN S :

z 7−→ µ(ρ)z + gb(z, z̄, ρ) (3.10)

iθ(ρ)
where λ(ρ) = (1 + ϕ(ρ))e
b with ϕb (ρNS ) = 0 and gb(z, z̄, ρ) is a smooth complex-valued function.
Applying Taylor expansion to the function gb, we obtain
X 1
gb(z, z̄, ρ) = gbkl (ρ)z k−l with gbkl ∈ C, k, l = 0, 1, . . . .
k! l!
k+l≥2
It is possible to define the Taylor coefficients by means of symmetric multi-linear vector functions.
gb20 (ρNS ) = hq2 , B(q1 , q1 )i,
gb11 (ρNS ) = hq2 , B(q1 , q¯1 )i,
(3.11)
gb02 (ρNS ) = hq2 , B(q¯1 , q¯1 )i,
gb21 (ρNS ) = hq2 , C(q1 , q1 , q¯1 )i.
The sign of the first Lyapunov coefficient l2 (ρN S ) determines the direction of NS bifurcation and is
defined by    
2
2 2
l2 (ρNS ) = Re λ22gb21 − Re (1−2λ 1 )λ2
2(1−λ1 ) g g11 | − 14 |b
b20 gb11 − 12 |b g02 | . (3.12)
According to the above discussion, the direction and stability of NS bifurcation can be presented in the
following theorem.

Theorem 3.2. Suppose (3.6) holds and l2 (ρN S ) 6= 0, then NS bifurcation at fixed point E(x∗ , y ∗ )
for system (1.10) if the ρ changes its value in small neighbourhood of N SE . Moreover, if l2 (ρN S ) <
0 (resp. l2 (ρN S ) > 0), then there is a smooth closed invariant curve that can bifurcate from E(x∗ , y ∗ ),
and the bifurcation is sub-critical (resp. super-critical).

4. Numerical Simulations
In this part, we will carry out numerical simulations to corroborate our theoretical findings for system
(1.10) that include bifurcation diagrams, phase portraits, MLEs and FDs.

Example 1: These are the chosen parameter values: a = 1.5, b = 0.5, α = 0.58 and, ρ fluctuates in
the range 0.3 ≤ ρ ≤ 0.5. Also the initial condition is (x0 , y0 ) = (1.998, 0.121). We identify a fixed point
E(x∗ , y ∗ ) = (2, 0.125) and bifurcation point for the system (1.10) is evaluated at ρ− = 0.349411.
and the eigenvalues of A(ρ− ) are λ1,2 = −1, 0.256747.
Let q1 , q2 ∈ R2 be two eigenvectors of A(ρ− ) and AT (ρ− ) corresponding to λ1,2 , respectively. There-
fore,
q1 ∼ (−1.43845, 1)T and q2 ∼ (0.179806, 1)T .
For hq1 , q2 i = 1 to told, we take normalized vector as q2 = γF q2 where, γF = 1.34887. Then we get
 CHAOTIC DYNAMICS OF THE FRACTIONAL ORDER SCHNAKENBERG MODEL AND ITS CONTROL 49

(a) (b)

(c) (d)

Figure 1. Flip Bifurcation diagram in (a) (ρ, x) plane, (b) (ρ, y) plane, (c) MLEs, (d) FDs .

q1 ∼ (−1.43845, 1)T and q2 ∼ (0.24253, 1.34887)T .

From (3.4), we obtain the Lyapunov coefficient l1 (ρ− ) = 3.93558 > 0. Therefore, the Flip bifurcation
is sub-critical and the requirements of Theorem 3.1 are established.

The bifurcation diagrams in Fig. 1 (a-b) demonstrate that fixed point E stability takes place for
ρ < ρ− , loses stability at ρ = ρ− , and a period doubling phenomenon results in chaos for ρ > ρ− . Fig.1
(c-d) displays the calculated MLEs and FDs associated to Fig.1(a-b).We see that different choices of
ρ result in the chaotic set with period −2, −4, −8, and −16 orbits. Dynamical states that are stable,
periodic, or chaotic are compatible with the sign in Fig.1(c-d), in accordance with the highest Lyapunov
exponent. For various values of ρ, the phase portraits of the bifurcation diagrams in Fig.1(a-b) are shown
in Fig.2.
Example 2: The following parameter values are chosen: a = 0.5, b = 0.5, α = 0.58 and, ρ fluctuates
in the range 0.45 ≤ ρ ≤ 0.6. Also the initial condition is (x0 , y0 ) = (0.998, 0.498). We obtain a fixed
point E(x∗ , y ∗ ) = (1, 0.5) and bifurcation point of the system (1.10) is evaluated at ρN S = 0.820228.
and the eigenvalues are λ1,2 = 0.5 ± 0.866025i. Also

d |λi (ρ)| A2e

|ρ=ρN S = − = 0.5 6= 0,
dρ 2
50 MD. JASIM UDDIN AND S. M. SOHEL RANA

(a) (b) (c)

(d) (e) (f)

(g) (h) (i)

Figure 2. Phase picture for various ρ values matching to Figure 1 a,b. Red ∗ is the fixed
point E0 .

A22e
−(T r(JE ))|ρ=ρN S 6= 0 ⇒ = 1 6= 2, 3.
A1e
Let q1 , q2 ∈ C2 be two complex eigenvectors corresponding to λ1,2 , respectively. Therefore,
q1 ∼ (−0.5 − 0.866025i, 1)T and q2 ∼ (0.5 − 0.866025i, 1)T .
For hq1 , q2 i = 1 to hold, q2 is normalized vector q2 = γN S q2 where, γN S = 0.5 − 0.288675i. Also, by
(3.11), the Taylor coefficients are
gb20 = 2.0 + 0.57735i, gb11 = 0.5 − 0.288675i, gb02 = 0.5 − 2.02073i, gb21 = −2 + 0.i.
From (3.12),The Lyapunov coefficient is obtained as l2 (ρN S ) = −0.75 < 0. As a result, the conditions
of Theorem 3.2 are met and the NS bifurcation is super-critical.

The NS bifurcation diagrams are shown in Fig.3(a,b), which reveals that the fixed point E is stable
while ρ < ρN S , but loses stability when ρ = ρN S , and exhibits an attractive closed invariant curve when
ρ > ρN S . Because of the presence of MLEs (3(c)), system dynamics are not stable.
The behavior of the smooth invariant curve as it separates from the stable fixed point and expands
in radius is nicely illustrated by the phase portraits (Fig.4) of the bifurcation diagrams in Fig.4 for
various values of rho.The closed curve abruptly vanishes as ρ increases, and for different values of ρ,
orbits with periods of −7, −10, −20, and −40 arise.
Example 3: As other parameter values change (e.g. parameter a), the Schnakenberg model may
exhibit more dynamic behavior in the Neimark-Sacker bifurcation diagram. A new Neimark-Sacker
bifurcation diagram is created when the parameter values are set as in Example 2 with ρ = 0.53145 and
a range between 0.1 ≤ a ≤ 0.3, as shown in Figure (5) (a-b). The system experiences a Neimark-Sacker
bifurcation at a = aN S = 0.11. The maximal Lyapunov exponent, which corresponds to Figure5(a-b),
 CHAOTIC DYNAMICS OF THE FRACTIONAL ORDER SCHNAKENBERG MODEL AND ITS CONTROL 51

(a) (b)

(c) (d)

Figure 3. NS Bifurcation diagram in (a) (ρ, x) plane, (b) (ρ, y) plane, (c) MLEs, (d) FDs.

is calculated and shown in Figure5(c), confirming the existence of chaos and the period window as a act
as a variable parameter. A thorough explanation of the behavior of the smooth invariant curve may be
found in the phase portrait of the bifurcation diagrams for different values of a illustrated in Figure6.
This diagram illustrates how the stable fixed point breaks off the smooth invariant curve as its radius
rises. Figure(7) displays the plot of the maximum Lyapunov exponents for two control parameters as a
2D projection onto the (ρ, a) and (ρ, b) plane. We note that the dynamics of system (1.10) shift from
chaotic to non-chaotic phase when the values of the control parameters a and b grow.
4.1. Fractal Dimension. The measurement of the fractal dimensions (FD) serves to identify a system’s
chaotic attractors and is defined by [10]
Pk
j=1 tj
D
bL = k + (4.1)
|tk+1 |
where k is the largest integer such that
k
X k+1
X
tj ≥ 0 and tj < 0,
j=1 j=1

and tj ’s are Lyapunov exponents.

Now, the fractal dimensions for the system (1.10) take the following form:
b L = 2 + t1 .
D (4.2)
|t2 |
52 MD. JASIM UDDIN AND S. M. SOHEL RANA

(a) (b) (c)

(d) (e) (f)

(g) (h) (i)

(j) (k)

Figure 4. Phase picture for various ρ values matching to Figure 3 a,b. Red ∗ is the fixed
point E0 .

Because the chaotic dynamics of the system (1.10) ( see Fig. 2) are quantified with the sign of
FD (see Fig. 3 (d)), it is certain that as the value of the parameter ρ increases, the dynamics of the
Fractional Order Schnakenberg system become unstable.

5. Chaos Control
It is believed that dynamical systems will be optimized in relation to some performance criterion and
will prevent chaos. In physics, biology, ecology, chemical engineering, telecommunications, and other
fields, chaotic behavior is studied. Additionally, the practical chaos management techniques can be
used in a variety of fields, including physics labs, biochemistry, turbulence, and cardiology, as well as
in communication systems. Recently, a lot of scholars have become quite interested in the problem of
managing chaos in discrete-time systems.
Controlling chaos is a challenging issue. For controlling chaos in fractional order Schnakenberg model
we introduce OGY [33], Hybrid control strategy [41] and state feedback[27] . In OGY method we can
not use ρ as a control parameter. We use a as a control parameter to implement OGY method.
 CHAOTIC DYNAMICS OF THE FRACTIONAL ORDER SCHNAKENBERG MODEL AND ITS CONTROL 53

(a) (b)

(c) (d)

Figure 5. NS Bifurcation diagram in (a) (a, x) plane, (b) (a, y) plane, (c) MLEs , (d) FDs

We can modify the system (1.10) as follows to apply OGY method:

ρα
a − xn + x2n yn = f˜1 (x, y, a),


xn+1 = xn +
Γ(α + 1)
(5.1)
ρα
b − x2n yn = f˜2 (x, y, a).


yn+1 = yn +
Γ(α + 1)

where a represents the chaos control parameter. In addition, it is presumable that |a − a0 | < ν̃, where
ν̃ > 0 and a0 represents the nominal parameter, lies in the chaotic areas. We employ a stabilizing
feedback control strategy to steer the trajectory toward the desired orbit. If the system (1.10) has
an unstable fixed point at (x+ , y + ) in a chaotic zone created by the formation of a Neimark-Sacker
bifurcation, the system (5.1) can be approximated in the area surrounding the unstable fixed point at
(x+ , y + ) by the following linear map:

xn+1 − x+ xn − x+
   
≈ A˜ee + B˜ee [a − a0 ] (5.2)
yn+1 − y + yn − y +
where
∂ f˜1 (x,y,a) ∂ f˜1 (x,y,a)
"  
ρα α
" # (−a+b)
#
ρ
1+ (a + b)2 Γ(α+1)
A˜ee = ∂x
∂ f˜2 (x,y,a)
∂y
∂ f˜2 (x,y,a)
= −2b
(a+b) Γ(α+1)
ρα ρ α .
∂x ∂y a+b Γ(α+1) 1 − (a + b)2 Γ(α+1)
54 MD. JASIM UDDIN AND S. M. SOHEL RANA

(a) (b) (c)

(d) (e) (f)

(g)

Figure 6. Phase picture for various a values matching to Figure 5 a,b. Red ∗ is
the fixed point E0 .

(a) (b)

Figure 7. Maximum Lyapunov exponents projected in two dimensions onto the (a) (ρ, a)
plane (b) (ρ, b) plane

and
∂ f˜1 (x,y,a) ρα
" # " #
B˜ee = ∂a
∂ f˜2 (x,y,a)
= Γ(α+1) .
∂a
0
Following that, we define the system’s 5.1 controllability matrix as follows:
 
ρα ρα ρα
 
(−a+b)
h i Γ(α+1) 1 + (a+b) Γ(α+1) Γ(α+1)
C˜ee = B˜ee : A˜ee B˜ee =   α 2 .
−2b ρ
0 a+b Γ(α+1)
 CHAOTIC DYNAMICS OF THE FRACTIONAL ORDER SCHNAKENBERG MODEL AND ITS CONTROL 55

The rank of C˜ee is then obvious to perceive to be 2. Considering that

xn − x+
 
[a − a0 ] = −K˜ee
yn − y +

where K˜ee = [σ˜e1 σ˜e2 ], system (5.1) becomes

xn+1 − x+ xn − x+
   
≈ [A˜ee − B˜ee K˜ee ] .
yn+1 − y + yn − y +
Additionally, the suitable controlled system of (1.10) is offered by

ρα
(a0 − σ˜e1 (xn − x+ ) − σ˜e2 (yn − y + )) − xn + x2n yn ,


xn+1 = xn +
Γ(α + 1)
(5.3)
ρα
b − x2n yn .


yn+1 = yn +
Γ(α + 1)
The fixed point (x+ , y + ) is additionally locally asymptotically stable if and only if both of the matrix’s
eigenvalues (A˜ee − B˜ee K˜ee ) are located within an open unit disk. Also,
"  
ρα ρα α
ρα
(−a+b)
#
ρ
1+ − Γ(α+1) σ˜
e1 (a + b)2 Γ(α+1) − Γ(α+1) σ˜
e2
A˜ee − B˜ee K˜ee = −2b
(a+b) Γ(α+1)
ρα ρ α .
a+b Γ(α+1) 1 − (a + b)2 Γ(α+1)

The Jacobian matrix (A˜ee − B˜ee K˜ee ) has the following characteristic equation:

ρα
   
(−a + b)
λe 2 − 2 − (a + b)2 + − σ˜e1 λe
(a + b) Γ(α + 1)
2 !
ρα ρα

1 3 2 2

3 3


+ a + b − a + a − b − 3a b + 3ab − b − (a + b) (5.4)
a+b Γ(α + 1) Γ(α + 1)
2 !
ρα α α
  
1 ρ ρ
+ (a + b) −1 + (a + b)2 σ˜e1 − 2bσ˜e2 = 0.
a+b Γ(α + 1) Γ(α + 1) Γ(α + 1)

If we consider eigenvalues λe1 and λe2 of the characteristic equation (5.4), we obtain

ρα
   
(−a + b)
λe1 + λe2 = 2 − (a + b)2 + − σ˜e1 ,
(a + b) Γ(α + 1)
ρα
 
1
a + b − a + a3 − b − 3a2 b + 3ab2


λe1 λe2 =
a+b Γ(α + 1)
!
2
ρα ρα ρα
   
1 3 3

2
+ − b − (a + b) + (a + b) −1 + (a + b) σ˜e1
a+b Γ(α + 1) Γ(α + 1) Γ(α + 1)
2 !
ρα

1
+ −2bσ˜e2 .
a+b Γ(α + 1)
(5.5)
Then, we must work out the solutions to the equations λe1 = ±1 and λe1 λe2 = 1 to get the lines of
marginal stability. Furthermore, these restrictions ensure that λe1 and λe2 are located inside the open
unit disk. Assume λe1 λe2 = 1 and from (5.5), we get
56 MD. JASIM UDDIN AND S. M. SOHEL RANA

ρα ρα
 
1 3 3
Ke1 = −(a + b) − (a + b) + (a + b)
a + b Γ(α + 1) Γ(α + 1)
2 !
ρα ρα ρα ρα
    
1 1
+ (a + b) −1 + (a + b)2 σ˜e1 − 2bσ˜e2 .
a + b Γ(α + 1) Γ(α + 1) a + b Γ(α + 1) Γ(α + 1)
Next, consider λe1 = 1, we obtain
2 !
ρα ρα

1 3 3
Ke2 = 4(a + b) − 2(a − b + (a + b) ) + (a + b)
a+b Γ(α + 1) Γ(α + 1)
2 !
ρα ρα ρα
  
1 2
+ (a + b) −2 + (a + b) σ˜e1 − 2bσ˜e2 .
a+b Γ(α + 1) Γ(α + 1) Γ(α + 1)
Lastly, if λe1 = −1 , then


ρα (a + b)3 + (a + b)3 σ˜e1 − 2bσ˜e2
Ke3 = − .
Γ(α + 1) a+b
Stable eigenvalues are then found in the triangle in the σ˜e1 , σ˜e2 plane enclosed by the straight lines
Ke1 , Ke2 , Ke3 for a given parametric value.
Hybrid control strategy is applied to system (1.10) to control chaos. We rewrite our uncontrolled
system (1.10) as
Xn+1 = G(Xn , δ) (5.6)
where δ ∈ R is bifurcation parameter and Xn ∈ R2 , G(.) is non-linear vector function. Applying hybrid
control strategy, the controlled system of (5.2) becomes

Xn+1 = ωe G(Xn , δ) + (1 − ωe )Xn (5.7)

where the control parameter is 0 < ωe < 1. Now, if we implement the above mentioned control strategy
to system (1.10), then we get the following controlled system

pα
 
a − xn + x2n yn + (1 − ωe )xn ,


xn+1 = ωe xn +
Γ(α + 1)
(5.8)
sα
 
2


yn+1 = ωe yn + b − xn yn + (1 − ωe )yn .
Γ(α + 1)
An approach called state feedback control is used to stabilize chaos at the stage of the system’s
(1.10) unstable trajectories. System (1.10) may be made to take on a controlled form by introducing a
feedback control law as the control force uee and is given below.

ρα
a − xn + x2n yn + uee ,


xn+1 = xn +
Γ(α + 1)
ρα (5.9)
b − x2n yn ,


yn+1 = yn +
Γ(α + 1)
uee = −k1 (xn − x+ ) − k2 (yn − y + ),
where (x+ , y + ) represents the positive fixed point of the system (1.10) and k1 and k2 stand for the
feedback gains.
Example 4: We use (a0 , b, α, ρ) = (0.15, 1.8, 0.58, 0.53145) to implement the OGY feedback control
technique for system (1.10).The unstable system (1.10) in this case has a single positive fixed point
 CHAOTIC DYNAMICS OF THE FRACTIONAL ORDER SCHNAKENBERG MODEL AND ITS CONTROL 57

(x+ , y + ) = (1.95, 0.473373). Then, in accordance with these parametric values, we provide the following
controlled system

xn+1 = xn + 0.777474 (0.15 − σ˜e1 (xn − 1.95) − σ˜e2 (yn − 0.473373)) − xn + x2n yn ,


(5.10)
yn+1 = yn + 0.777474 1.8 − x2n yn .

where K̃ = [σ˜e1 σ˜e2 ] serve as the gain matrix. We also get,

 
˜ 1.65786 2.95634
Aee = ,
−1.43534 −1.95634
 
0.777474
B˜ee = ,
0
 
0.777474 1.28895
C˜ee = .
0 −1.11594
The rank of the matrix C˜ee is therefore easily verified to be 2. Consequently, it is possible to control
the system (5.10) and the jacobian matrix of the controlled system is given by
 
1.65786 − 0.777474σ˜e1 2.95634 − 0.777474σ˜e2
A˜ee − B˜ee K˜ee = .
−1.43534 −1.95634
For marginal stability, the lines Ke1 , Ke2 and Ke3 are provided by
Ke1 = −0.000000986495 + 1.52101σ˜e1 − 1.11594σ˜e2 = 0,
Ke2 = 1.70152 + 0.743533σ˜e1 − 1.11594σ˜e2 = 0,
Ke3 = −2.29848 − 2.29848σ˜e1 + 1.11594σ˜e2 = 0.
The marginal lines Ke1 , Ke2 and Ke3 define the stable triangular area for the controlled system (5.10)
is displayed in Figure(8).
Determining the hybrid control technique’s efficiency in reducing chaotic (unstable) system dynamics,
we are taken the parameter values discussed for OGY method except ρ = 0.5812 > ρN S . As a result,
it demonstrates that the fixed point E(1.95, 0.473373) of system (1.10) is unstable, however for the
controlled system(5.8), this fixed point is stable iff 0 < ωe < 0.9494221636969571. By taking ωe = 0.85,
which demonstrates that the fixed point E is a sink for the controlled system(5.8) and removes the
unstable system dynamics around E. The stable region and stable trajectories are shown in Figure(8).
We have carried out numerical simulations (in Figure(8)) to examine the operation of the feedback
control method’s control of chaos in an unstable condition. The parameters values will be same as we
choose for OGY method. The feedback gains are choosen as k1 = −0.42 and k2 = −0.35.

6. Conclusions
In this study, a novel fractional order Schnakenberg model is discussed. The Caputo fractional
derivative idea is the basis for such a fractional order model. We show that the system (1.10) can
experience a bifurcation (flip or NS) at a specific positive fixed point E if fluctuations about the sets
P DE or N SE . The center manifold theorem and bifurcation theory are used to achieve this. The model
displays a range of complex dynamical behaviors as ρ and α are changed, including the appearance of
flip and NS bifurcations, period 2, 4, 6, 7, 8, 10, 16, 20 and 40 orbits, and quasi-periodic orbits, as well as
attracting invariant circles and chaotic sets. By computing the maximal Lyapunov exponents and fractal
dimension, we can verify the existence of chaos. Typically, one of the traits that suggests the existence
of chaos is a positive Lyapunov exponent.We demonstrate that the values of the Lyapunov exponent
58 MD. JASIM UDDIN AND S. M. SOHEL RANA

(i) (ii) (iii)

(iv) (v) (vi)

Figure 8. (i-iii) Stable region in OGY method, Hybrid method and state feedback
methods, (iv-vi) Stable system trajectories

can alternately be negative and positive. By creating pathways from periodic and quasi-periodic states,
the two bifurcations lead the system to quickly transition from steady state to chaotic dynamical
behavior. Alternatively, chaotic dynamics might occur or vanish concurrently with the appearance of
bifurcations. Finally, we reduce unstable system trajectories by employing an OGY, Hybrid, and state
feedback control technique. Though it remains a challenging topic, investigating multiple parameter
bifurcation in the system.

References
[1] M. A. M. Abdelaziz, A. I. Ismail, F. A. Abdullah and M. H. Mohd, Bifurcations and chaos in a discrete SI epidemic
model with fractional order, Advances in Difference Equations, 2018(2018), 1-19.
[2] R. P. Agarwal and P. J. Y. Wong, Advanced Topics in Difference Equations, Mathematics and its Applications 404,
Kluwer Academic Publishers Group, Dordrecht, 1997.
[3] R. P. Agarwal, A. M. A. El-Sayed and S. M. Salman, Fractional-order Chua’s system: discretization, bifurcation and
chaos, Advances in Difference Equations, 1(2013), 1–13.
[4] W. M. Ahmad and J. C. Sprott, Chaos in fractional-order autonomous nonlinear systems, Chaos, Solitons & Fractals,
16(2003), 339–351.
[5] K. S. Al Noufaey, Semi-analytical solutions of the Schnakenberg model of a reaction-diffusion cell with feedback,
Results in Physics, 9(2018), 609-614.
[6] I. Ameen and P. Novati, The solution of fractional order epidemic model by implicit Adams methods, Applied
Mathematical Modelling, 43(2017),78–84.
[7] R. L. Bagley and R. A. Calico, Fractional order state equations for the control of viscoelasticallydamped structures,
Journal of Guidance, Control, and Dynamics, 14(1991), 304–311.
[8] E. Balcı, İ. Öztürk and S. Kartal , Dynamical behaviour of fractional order tumor model with Caputo and conformable
fractional derivative, Chaos, Solitons & Fractals, 123(2019), 43–51.
[9] E. Balci, S. Kartal and I. Ozturk , Comparison of dynamical behavior between fractional order delayed and discrete
conformable fractional order tumor-immune system, Mathematical Modelling of Natural Phenomena, 16(2021):3.
[10] J. H. E. Cartwright, Nonlinear stiffness, Lyapunov exponents, and attractor dimension, Physics Letters A, 264(1999),
298-302.
[11] Q. Din, A novel chaos control strategy for discrete-time Brusselator models, Journal of Mathematical Chemistry,
56(2018), 3045–3075.
[12] Z. F. El Raheem and S. M. Salman, On a discretization process of fractional-order logistic differential equation,
Journal of the Egyptian Mathematical Society, 22(2014), 407–412.
 CHAOTIC DYNAMICS OF THE FRACTIONAL ORDER SCHNAKENBERG MODEL AND ITS CONTROL 59

[13] A. A. Elsadany and A. E. Matouk, Dynamical behaviors of fractional-order Lotka–Volterra predator–prey model and
its discretization, Journal of Applied Mathematics and Computing, 49(2015), 269–283.
[14] C. A. Floudas and X. Lin, Continuous-time versus discrete-time approaches for scheduling of chemical processes: a
review, Computers & Chemical Engineering, 28(2004), 2109–2129.
[15] Z. Hammouch and T. Mekkaoui, Numerical simulations for a variable order fractional Schnakenberg model, AIP
Conference Proceedings, 1637(2014), 1450–1455.
[16] C. Huang, J. Cao, M. Xiao, A. Alsaedi and T. Hayat, Bifurcations in a delayed fractional complex-valued neural
network, Applied Mathematics and Computation, 292(2017), 210–227.
[17] C. Huang, J. Cao, M. Xiao and A. Alsaedi, Controlling bifurcation in a delayed fractional predator–prey system with
incommensurate orders, Applied Mathematics and Computation, 293(2017), 293–310.
[18] C. Huang, J. Cao, M. Xiao, A. Alsaedi and T. Hayat, Effects of time delays on stability and Hopf bifurcation in
a fractional ring-structured network with arbitrary neurons, Communications in Nonlinear Science and Numerical
Simulation,57 (2018), 1–13.
[19] C. Huang, Y. Meng, J. Cao and A. Alsaedi, New bifurcation results for fractional BAM neural network with leakage
delay, Chaos, Solitons & Fractals, 100(2017), 31–44.
[20] M. Ichise, Y. Nagayanagi and T. Kojima,An analog simulation of non-integer order transfer functions for analysis
of electrode processes, Journal of Electroanalytical Chemistry and Interfacial Electrochemistry, 33(1971), 253–265.
[21] H. Jafari and V. Daftardar-Gejji, Solving a system of nonlinear fractional differential equations using Adomian
decomposition, Journal of Computational and Applied Mathematics, 196(2006),644–651.
[22] R. Kapral , Discrete models for chemically reacting systems, Journal of Mathematical Chemistry, 6(1991), 113–163.
[23] F. M. Khan, A. Ali, N. Hamadneh, Abdullah, M. N. Alam , Numerical Investigation of Chemical Schnakenberg
Mathematical Model, Journal of Nanomaterials, 2021(2021): Article ID 9152972.
[24] A. Q. Khan and T. Khalique, Neimark-Sacker bifurcation and hybrid control in a discrete-time Lotka-Volterra model,
Mathematical Methods in the Applied Sciences, 43(2020), 5887–5904.
[25] A. Q. Khan, S. A. H. Bukhari and M.B. Almatrafi, Global dynamics,Neimark-Sacker bifurcation and hybrid control
in a Leslie’s prey-predator model, Alexandria Engineering Journal, 61(2022):11391–11404.
[26] Y. A. Kuznetsov, I. A. Kuznetsov and Y. Kuznetsov, Elements of Applied Bifurcation Theory, Springer-Verlag, New
York, 1998.
[27] S. Lynch, Dynamical Systems with Applications Using Mathematica, 2nd ed, Springer International Publishing AG,
2007.
[28] P. Liu, J. Shi, Y. Wang and X. Feng, Bifurcation analysis of reaction–diffusion Schnakenberg model,Journal of
Mathematical Chemistry, 51(2013), 2001–2019.
[29] S.A. Manaa, Some numerical methods for Schnackenberg model, International Journal of Engineering Inventions,
2(2013), 71–78.
[30] R. Magin, M.D. Ortigueira, I. Podlubny and J. Trujillo, On the fractional signals and systems,Signal Processing,
91(2011), 350-371.
[31] S. Momani and Z. Odibat, Numerical approach to differential equations of fractional order, Journal of Computational
and Applied Mathematics, 207(2007),96–110.
[32] J. D. Murray, Mathematical Biology: I. An Introduction, Springer, New York, 2002.
[33] E. Ott, C. Grebogi and J. A. Yorke, Controlling chaos, Physical Review Letters, 64(1990),1196–1199.
[34] R. K. Pearson , Discrete-time Dynamic Models, Oxford University Press, Oxford, 1999.
[35] I. Podlubny, Fractional Differential Equations, Academic Press, New York,1999.
[36] S. M. Rana and U. Kulsum, Bifurcation analysis and chaos control in a discrete-time predator-prey system of Leslie
type with simplified Holling type IV functional response, Discrete Dynamics in Nature and Society, 2017(2017):
Article ID 9705985.
[37] S. M. Rana, Dynamics and chaos control in a discrete-time ratio-dependent Holling-Tanner model, Journal of the
Egyptian Mathematical Society, 27(2019),1–16.
[38] J. Schnakenberg, Simple chemical reaction systems with limit cycle behavior, Journal of Theoretical Biology,
81(1979), 389-400.
[39] G. Wen, Criterion to identify Hopf bifurcations in maps of arbitrary dimension, Physical Review E, 72(2005):
026201.
[40] S. Yao , New bifurcation critical criterion of Flip-Neimark-Sacker bifurcations for two-parameterized family of n-
dimensional discrete systems, Discrete Dynamics in Nature and Society, 2012(2012): Article ID 264526.
[41] L. G. Yuan and Q. G. Yang, Bifurcation, invariant curve and hybrid control in a discrete-time predator–prey system,
Applied Mathematical Modelling, 39(2015), 2345–2362.
60 MD. JASIM UDDIN AND S. M. SOHEL RANA

Md. Jasim Uddin, Corresponding author, Department of Mathematics, University of Dhaka, Dhaka-1000,
Bangladesh.
Email address: jasim.uddin@du.ac.bd

S. M. Sohel Rana, Department of Mathematics, University of Dhaka, Dhaka-1000, Bangladesh.

Email address: srana.mthdu@gmail.com