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 https://ojs.lib.uwo.ca/mase https://dx.doi.org/https://doi.org/10.5206/mase/15355 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: X m1 m−1 A,B m2→ Y, 2X + Y m3→ 3X, 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 dτ = m1NA −m−1NX + m3N2XNY , dNY dτ = m2NB −m3N2XNY , (1.1) 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 + x2y, ẏ = b−x2y, (1.2) where x = √ m3 m−1 NX, y = √ m3 m−1 NY , τ −→ t = τm−1, a = m1 m−1 √ m3 m−1 NA, b = √ m2 m−1 NB. 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), Dαy(t) = b−x2(t)y(t), (1.3) 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: Dαx(t) = a−x([ t ρ ]ρ) + x2([ t ρ ]ρ)y([ t ρ ]ρ), Dαy(t) = b−x2([ t ρ ]ρ)y([ t ρ ]ρ). (1.4) First, let t ∈ [0,ρ), so t ρ ∈ [0, 1). Thus, we obtain Dαx(t) = a−x0 + x20y0; Dαy(t) = b−x20y0. (1.5) The solution of (1.5) is simplified to x1(t) = x0 + J α ( a−x0 + x20y0 ) = x0 + tα αΓ(α) ( a−x0 + x20y0 ) , y1(t) = y0 + J α ( b−x20y0 ) = y0 + tα αΓ(α) ( b−x20y0 ) . (1.6) Second, let t ∈ [ρ, 2ρ), so t ρ ∈ [1, 2). Then Dαx(t) = a−x1 + x21y1 Dαy(t) = b−x21y1 (1.7) which have the following solution x2(t) = x1(ρ) + J α ρ ( a−x1 + x21y1 ) = x1(ρ) + (t−ρ)α αΓ(α) ( a−x1 + x21y1 ) , y2(t) = y1(ρ) + J α ρ ( b−x21y1 ) = y1(ρ) + (t−ρ)α αΓ(α) ( b−x21y1 ) , (1.8) where Jαρ ≡ 1 Γ(α) ∫ t ρ (t− τ)α−1dτ, α > 0. The result of doing the discretization process n times over xn+1(t) = xn(nρ) + (t−nρ)α αΓ(α) ( a−xn(nρ) + x2n(nρ)yn(nρ) ) , yn+1(t) = yn(nρ) + (t−nρ)α αΓ(α) ( b−x2n(nρ)yn(nρ) ) , (1.9) 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 xn+1 = xn + ρα Γ(α + 1) ( a−xn + x2nyn ) , yn+1 = yn + ρα Γ(α + 1) ( b−x2nyn ) . (1.10) 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: J (x∗,y∗) =   ( 1 + (−1 + 2x∗y∗) ρ α Γ(α+1) ) −x∗ 2 ρα Γ(α+1) −2x∗y∗ ρ α Γ(α+1) ( 1 −x∗ 2 )   . (2.1) Now at E(x∗,y∗), the Jacobian matrix is given by JE = ( 1 + (−a+b) a+b ρα Γ(α+1) (a + b)2 ρ α Γ(α+1) − 2b a+b ρα Γ(α+1) ( 1 − (a + b)2 ρ α Γ(α+1) ) ) . (2.2) The characteristic polynomial of the Jacobian Matrix can be written as F̃(λ) := λ2 − Tr(JE) λ + Det(JE) = 0, (2.3) where Tr(JE) is the trace and Det(JE) is the determinant of the Jacobian matrix JE, and is given by Tr(JE) =2 + (−a + b) a + b − (a + b)2 ρα Γ(α + 1) , Det(JE) = a3(−1 + ρ̂)ρ̂ + 3a2b(−1 + ρ̂)ρ̂ + a(−1 + ρ̂)(−1 + 3b2ρ̂) + b(1 + ρ̂− b2ρ̂ + b2ρ̂2) a + b . (2.4) where ρ̂ = ρα/Γ(α + 1). 44 MD. JASIM UDDIN AND S. M. SOHEL RANA The eigenvalues of the system(2.3) can be written as λ1,2 = Tr(JE) ± √ (Tr(JE))2 − 4Det(JE) 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 PDE =  (a,b,ρ,α) : ρ = ( Γ(1 + α) −A2e ± √ L A1e ) 1 α = ρ±,L ≥ 0   where A1e = (a + b) 2, A2e = − a + a3 − b + 3a2b + 3ab2 + b3 a + b = − a− b + (a + b)3 a + b = − A4e a + b , A3e = 4 A4e = a− b + (a + b)3, L = A22e −A1eA3e. The system (1.10), however, undergoes a flip bifurcation at E when the parameters (a,b,ρ,α) fluctuate within a narrow region of PDE. Also let NSE = { (a,b,ρ,α) : ρ = ( Γ(1 + α) −A2e A1e ) 1 α = ρNS,L < 0 } . If the parameters (a,b,ρ,α) vary around the set NSE, 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, ρ < ρNS(stable focus), − source (i)L ≥ 0, ρ > ρ+(unstable node), (ii)L < 0, ρ > ρNS(unstable focus), − non-hyperbolic (i)L ≥ 0, ρ = ρ−or ρ = ρ+( saddle with flip), (ii)(i)L < 0, ρ = ρNS( 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 PDE. Consider the system (1.10) at equilibrium point E(x∗,y∗). Assume the parameters lie in PDE. Let L ≥ 0 and ρ = ρ− = ( Γ(1 + α) −A2e − √ L A1e ) 1 α . 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̂ ŷ ) → A(ρ−) ( x̂ ŷ ) + ( F1(x̂, ŷ,ρ−) F2(x̂, ŷ,ρ−) ) (3.2) where X = (x̂, ŷ)T and F1(x̂, ŷ,ρ−) = 1 (a + b)5 [ a− b + (a + b)3 − √ L ] x̂ ( 2b3ŷ + a2(2a + x̂)ŷ + b2(6a + x̂)ŷ ) + 1 (a + b)5 [ a− b + (a + b)3 − √ L ] x̂ ( b(x̂ + 6a2ŷ + 2ax̂ŷ) ) , F2(x̂, ŷ,ρ−) = − 1 (a + b)5 [ a− b + (a + b)3 − √ L ] x̂ ( 2b3ŷ + a2(2a + x̂)ŷ + b2(6a + x̂)ŷ ) − 1 (a + b)5 [ a− b + (a + b)3 − √ L ] x̂ ( b(x̂ + 6a2ŷ + 2ax̂ŷ) ) . (3.3) The system(3.2) can be expressed as Xn+1 = AXn + 1 2 B (Xn,Xn) + 1 6 C (Xn,Xn,Xn) + O ( ‖Xn‖ 4 ) where B(x,y) = ( B1(x,y) B2(x,y) ) and C(x,y,u) = ( C1(x,y,u) C2(x,y,u) ) are symmetric multi-linear vector functions of x,y,u ∈ R2 and defined as follows: B1(x,y) = 2∑ j,k=1 δ2F1(ξ,ρ) δξjδξk ∣∣∣∣∣∣ ξ=0 xjyk = 1 (a + b)5 [ a− b + (a + b)3 − √ L ]( a3(x2y1 + x1y2) ) + 1 (a + b)5 ( 3a2b(x2y1 + x1y2) ) + 1 (a + b)5 ( b3(x2y1 + x1y2) + 3a 2b(x2y1 + x1y2) + bx1y1 ) , B2(x,y) = 2∑ j,k=1 δ2F2(ξ,ρ) δξjδξk ∣∣∣∣∣∣ ξ=0 xjyk = 1 (a + b)5 [ a− b + (a + b)3 − √ L ]( (a + b)3(x2y1 + x1y2) + bx1y1 ) , 46 MD. JASIM UDDIN AND S. M. SOHEL RANA and C1(x,y,u) = 2∑ j,k,l=1 δ2F1(ξ,ρ) δξjδξkδξl ∣∣∣∣∣∣ ξ=0 xjykul = 1 (a + b)3 ( a− b + (a + b)3 − √ L ) (u1x1y2) + 1 (a + b)3 ( a− b + (a + b)3 − √ L ) (u1x1y2) , C2(x,y,u) = 2∑ j,k,l=1 δ2F1(ξ,ρ) δξjδξkδξl ∣∣∣∣∣∣ ξ=0 xjykul = − 1 (a + b)3 ( a− b + (a + b)3 − √ L ) (u2x1y1 + u1x2y1) − 1 (a + b)3 ( a− b + (a + b)3 − √ L ) (u1x1y2) . 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 q1 = ( −(a+b) 3(−a−3b−+(a+b)3− √ L) 2b(a−b+(a+b)3− √ L) 1 ) = ( q11 1 ) , q2 = ( −a+3b−(a+b) 3+ √ L a−b+(a+b)3− √ L 1 ) = ( q21 1 ) . In order to get 〈q1,q2〉 = 1, where 〈q1,q2〉 = q11q21 + q12q22, we have to use the normalized vector q2 = γFq2, with γF = 1/(1 + q11p11). 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 l1 (ρ−) = 1 6 〈q2,C(q1,q1,q1)〉− 1 2 〈 q2,B ( q2, (A− I)−1B(q1,q1) )〉 . (3.4) 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 PDE. 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,α,ρ) ∈ NSE and L < 0, then the eigenvalues of system (2.3) are given by λ,λ̄ = Tr(JE) ± i √ 4Det(JE) −Tr(JE)2 2 (3.5) where Tr(JE) = 2 − (a− b)(a− b + (a + b)3) (a + b)4 − a− b + (a + b)3 a + b , Det(JE) = 1. Let ρNS = −A2e/A1e. Then, the transversality and non-resonance conditions respectively read d|λi(ρ)| dρ |ρ=ρNS = − A2e 2 6= 0; −(Tr(JE))|ρ=ρNS 6= 0 ⇒ A22 A1e 6= 2, 3. (3.6) 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 F1(x̂, ŷ,ρNS) = (a− b + (a + b)3)x̂ ( 2b3ŷ + a2(2a + x̂)ŷ + b2(6a + x̂)ŷ + b(x̂ + 6a2ŷ + 2ax̂ŷ) ) (a + b)5 , F2(x̂, ŷ,ρNS) = − (a− b + (a + b)3)x̂ ( 2b3ŷ + a2(2a + x̂)ŷ + b2(6a + x̂)ŷ + b(x̂ + 6a2ŷ + 2ax̂ŷ) ) (a + b)5 . (3.8) The system(3.2) can be expressed as Xn+1 = AXn + 1 2 B (Xn,Xn) + 1 6 C (Xn,Xn,Xn) + O ( ‖Xn‖ 4 ) where B1(x,y) = 2(a− b + (a + b)3) (a + b)5 [ (a + b)3(x2y1 + x1y2) + bx1y1 ] , B2(x,y) = − 2(a− b + (a + b)3) (a + b)5 [ (a + b)3(x2y1 + x1y2) + bx1y1 ] , and C1(x,y,u) = 2(a− b + (a + b)3) (a + b)5 (u2x1y1 + u1x2y1 + u1x1y2) , C2(x,y,u) = − 2(a− b + (a + b)3) (a + b)5 (u2x1y1 + u1x2y1 + u1x1y2) . Let q1,q2 ∈ C2 be two eigenvectors of A(ρNS) and AT (ρNS)for eigenvalue λ(ρNS) and λ̄(ρNS) such that A (ρNS) q1 = λ (ρNS) q1, A (ρNS) q̄1 = λ̄ (ρNS) q̄1, AT (ρNS) q2 = λ̄ (ρNS) q2, A T (ρNS) q̄2 = λ (ρNS) q̄2. (3.9) Therefore, using a straight calculation, we obtain q1 = ( −−a+b+(a+b) 3+ √ L 4b 1 ) = ( ζ1 + iζ2 1 ) , q2 = ( −a+b+(a+b)3− √ L 2(a+b)3 1 ) = ( ξ1 + iξ2 1 ) . where ζ1 = − −a + b + (a + b)3 4b , ζ2 = − √ −L 4b , ξ1 = −a + b + (a + b)3 2(a + b)3 , ξ2 = − √ −L 2(a + b)3 . For 〈q1,q2〉 = 1 to hold, we can take normalized vector as q2 = γNSq2 where γNS = 1 1 + (ξ1 + iξ2)(ζ1 − iζ2) . The following is how the eigenvectors are calculated: 48 MD. JASIM UDDIN AND S. M. SOHEL RANA q1 = ( ζ1 + iζ2 1 ) , q2 = ( ξ1+iξ2 1+(ξ1+iξ2)(ζ1−iζ2) 1 1+(ξ1+iξ2)(ζ1−iζ2) ) . We break down X ∈ R2 as X = zq1 + z̄q̄1 by considering ρ vary near to ρNS and for z ∈ C. The precise formulation of z is z = 〈q2,X〉. As a result, the system (2.3) changed to the following system for |ρ| near ρNS: z 7−→ µ(ρ)z + ĝ(z, z̄,ρ) (3.10) where λ(ρ) = (1 + ϕ̂(ρ))eiθ(ρ) with ϕ̂ (ρNS) = 0 and ĝ(z, z̄,ρ) is a smooth complex-valued function. Applying Taylor expansion to the function ĝ, we obtain ĝ(z, z̄,ρ) = ∑ k+l≥2 1 k! l! ĝkl(ρ)z k−l with ĝkl ∈ C,k, l = 0, 1, . . . . It is possible to define the Taylor coefficients by means of symmetric multi-linear vector functions. ĝ20 (ρNS) = 〈q2,B(q1,q1)〉, ĝ11 (ρNS) = 〈q2,B(q1, q̄1)〉, ĝ02 (ρNS) = 〈q2,B(q̄1, q̄1)〉, ĝ21 (ρNS) = 〈q2,C(q1,q1, q̄1)〉. (3.11) The sign of the first Lyapunov coefficient l2(ρNS) determines the direction of NS bifurcation and is defined by l2 (ρNS) = Re ( λ2ĝ21 2 ) − Re ( (1−2λ1)λ22 2(1−λ1) ĝ20ĝ11 ) − 1 2 |ĝ11| 2 − 1 4 |ĝ02| 2 . (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(ρNS) 6= 0, then NS bifurcation at fixed point E(x∗,y∗) for system (1.10) if the ρ changes its value in small neighbourhood of NSE. Moreover, if l2(ρNS) < 0 (resp. l2(ρNS) > 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 〈q1,q2〉 = 1 to told, we take normalized vector as q2 = γFq2 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 ρNS = 0.820228. and the eigenvalues are λ1,2 = 0.5 ± 0.866025i. Also d |λi(ρ)| dρ |ρ=ρNS = − A2e 2 = 0.5 6= 0, 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. −(Tr(JE))|ρ=ρNS 6= 0 ⇒ A22e A1e = 1 6= 2, 3. 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 〈q1,q2〉 = 1 to hold, q2 is normalized vector q2 = γNSq2 where, γNS = 0.5 − 0.288675i. Also, by (3.11), the Taylor coefficients are ĝ20 = 2.0 + 0.57735i, ĝ11 = 0.5 − 0.288675i, ĝ02 = 0.5 − 2.02073i, ĝ21 = −2 + 0.i. From (3.12),The Lyapunov coefficient is obtained as l2(ρNS) = −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 ρ < ρNS, but loses stability when ρ = ρNS, and exhibits an attractive closed invariant curve when ρ > ρNS. 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 = aNS = 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] D̂L = k + ∑k j=1 tj |tk+1| (4.1) where k is the largest integer such that k∑ j=1 tj ≥ 0 and k+1∑ j=1 tj < 0, and tj’s are Lyapunov exponents. Now, the fractal dimensions for the system (1.10) take the following form: D̂L = 2 + t1 |t2| . (4.2) 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: xn+1 = xn + ρα Γ(α + 1) ( a−xn + x2nyn ) = f̃1(x,y,a), yn+1 = yn + ρα Γ(α + 1) ( b−x2nyn ) = f̃2(x,y,a). (5.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+ yn+1 −y+ ] ≈ Ãee [ xn −x+ yn −y+ ] + B̃ee [a−a0] (5.2) where Ãee = [ ∂f̃1(x,y,a) ∂x ∂f̃1(x,y,a) ∂y ∂f̃2(x,y,a) ∂x ∂f̃2(x,y,a) ∂y ] = [ ( 1 + (−a+b) (a+b) ρα Γ(α+1) ) (a + b)2 ρ α Γ(α+1) −2b 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 B̃ee = [ ∂f̃1(x,y,a) ∂a ∂f̃2(x,y,a) ∂a ] = [ ρα Γ(α+1) 0 ] . Following that, we define the system’s 5.1 controllability matrix as follows: C̃ee = [ B̃ee : ÃeeB̃ee ] =   ραΓ(α+1) ( 1 + (−a+b) (a+b) ρα Γ(α+1) ) ρα Γ(α+1) 0 −2b a+b ( ρα Γ(α+1) )2   . 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 [a−a0] = −K̃ee [ xn −x+ yn −y+ ] where K̃ee = [σ̃e1 σ̃e2], system (5.1) becomes [ xn+1 −x+ yn+1 −y+ ] ≈ [Ãee − B̃eeK̃ee] [ xn −x+ yn −y+ ] . Additionally, the suitable controlled system of (1.10) is offered by xn+1 = xn + ρα Γ(α + 1) ( (a0 − σ̃e1(xn −x+) − σ̃e2(yn −y+)) −xn + x2nyn ) , yn+1 = yn + ρα Γ(α + 1) ( b−x2nyn ) . (5.3) The fixed point (x+,y+) is additionally locally asymptotically stable if and only if both of the matrix’s eigenvalues (Ãee − B̃eeK̃ee) are located within an open unit disk. Also, Ãee − B̃eeK̃ee = [ ( 1 + (−a+b) (a+b) ρα Γ(α+1) − ρ α Γ(α+1) σ̃e1 ) (a + b)2 ρ α Γ(α+1) − ρ α Γ(α+1) σ̃e2 −2b a+b ρα Γ(α+1) 1 − (a + b)2 ρ α Γ(α+1) ] . The Jacobian matrix (Ãee − B̃eeK̃ee) has the following characteristic equation: λe 2 − ( 2 − ( (a + b)2 + (−a + b) (a + b) − σ̃e1 ) ρα Γ(α + 1) ) λe + 1 a + b ( a + b− ( a + a3 − b− 3a2b + 3ab2 ) ρα Γ(α + 1) − ( b3 − (a + b)3 )( ρα Γ(α + 1) )2) + 1 a + b ( (a + b) ρα Γ(α + 1) ( −1 + (a + b)2 ρα Γ(α + 1) ) σ̃e1 − 2bσ̃e2 ( ρα Γ(α + 1) )2) = 0. (5.4) If we consider eigenvalues λe1 and λe2 of the characteristic equation (5.4), we obtain λe1 + λe2 = ( 2 − ( (a + b)2 + (−a + b) (a + b) − σ̃e1 ) ρα Γ(α + 1) ) , λe1λe2 = 1 a + b ( a + b− ( a + a3 − b− 3a2b + 3ab2 ) ρα Γ(α + 1) ) + 1 a + b ( − ( b3 − (a + b)3 )( ρα Γ(α + 1) )2 + (a + b) ρα Γ(α + 1) ( −1 + (a + b)2 ρα Γ(α + 1) ) σ̃e1 ) + 1 a + b ( −2bσ̃e2 ( ρα Γ(α + 1) )2) . (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 Ke1 = 1 a + b ρα Γ(α + 1) ( −(a + b) − (a + b)3 + (a + b)3 ρα Γ(α + 1) ) + 1 a + b ρα Γ(α + 1) ( (a + b) ( −1 + (a + b)2 ρα Γ(α + 1) ) σ̃e1 ) − 1 a + b ρα Γ(α + 1) ( 2bσ̃e2 ( ρα Γ(α + 1) )2) . Next, consider λe1 = 1, we obtain Ke2 = 1 a + b ( 4(a + b) − 2(a− b + (a + b)3) ρα Γ(α + 1) + (a + b)3 ( ρα Γ(α + 1) )2) + 1 a + b ( (a + b) ρα Γ(α + 1) ( −2 + (a + b)2 ρα Γ(α + 1) ) σ̃e1 − 2bσ̃e2 ( ρα Γ(α + 1) )2) . Lastly, if λe1 = −1 , then Ke3 = − ρα Γ(α + 1) ( (a + b)3 + (a + b)3σ̃e1 − 2bσ̃e2 ) 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 = ωeG(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 xn+1 = ωe ( xn + pα Γ(α + 1) ( a−xn + x2nyn )) + (1 −ωe)xn, yn+1 = ωe ( yn + sα Γ(α + 1) ( b−x2nyn )) + (1 −ωe)yn. (5.8) 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. xn+1 = xn + ρα Γ(α + 1) ( a−xn + x2nyn ) + uee, yn+1 = yn + ρα Γ(α + 1) ( b−x2nyn ) , uee = −k1(xn −x+) −k2(yn −y+), (5.9) 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 + x2nyn ) , yn+1 = yn + 0.777474 ( 1.8 −x2nyn ) . (5.10) where K̃ = [σ̃e1 σ̃e2] serve as the gain matrix. We also get, Ãee = [ 1.65786 2.95634 −1.43534 −1.95634 ] , B̃ee = [ 0.777474 0 ] , C̃ee = [ 0.777474 1.28895 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 Ãee − B̃eeK̃ee = [ 1.65786 − 0.777474σ̃e1 2.95634 − 0.777474σ̃e2 −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 > ρNS. 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 PDE or NSE. 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 1. Introduction 2. Stability of Fixed point 3. Bifurcation Analysis 3.1. Flip Bifurcation 3.2. Neimark-Sacker Bifurcation 4. Numerical Simulations 4.1. Fractal Dimension. 5. Chaos Control 6. Conclusions References