Stability Results for Two-Dimensional Systems of Fractional-Order Difference Equations

Linear autonomous incommensurate systems that consist of two fractional-order difference equations of Caputo-type are studied in terms of their asymptotic stability and instability properties. More precisely, the asymptotic stability of the considered linear system is fully characterized, in terms of the fractional orders of the considered Caputo-type differences, as well as the elements of the linear system’s matrix and the discretization step size. Moreover, fractional-order-independent sufficient conditions are also derived for the instability of the system under investigation. With the aim of exemplifying the theoretical results, a fractional-order discrete version of the FitzHugh–Nagumo neuronal model is constructed and analyzed. Furthermore, numerical simulations are undertaken in order to substantiate the theoretical findings, showing that the membrane potential may exhibit complex bursting behavior for suitable choices of the model parameters and fractional orders of the Caputo-type differences.


Introduction
The domain of fractional calculus is regarded as an important tool in modelling various phenomena from different scientific and engineering fields [1][2][3][4][5][6]. Many recent research results conjecture that fractional-order systems lead to more precise results in a wide range of practical applications, relying on the fact that memory and hereditary properties of different processes may be successfully modeled while using fractional-order derivatives.
Stability and linearization results have been recently explored for continuous-time fractional-order systems in [7,8], as well as for discrete-time fractional systems [9,10]. The solutions of initial value problems that involve fractional-order Caputo-type and Riemann-Liouville-type difference equations with positive orders were given in [11], but a qualitative study regarding the stability properties of these types of equations or systems of several equations have not been explored. Stability results for linear systems with Caputo fractional-order difference equations with variable order of convolution type have been studied in [12], along with the recurrence formulas for the solutions of linear initial value problems for the considered fractional operators. Still, when compared to fractional-order differential equations [13][14][15][16][17][18], their discrete-time counterparts, namely fractional-order difference equations, have not received the same amount of attention, especially regarding their qualitative theory and, in particular, their stability and instability properties.
is called the forward h-difference operator.
Definition 3. Let q ∈ (0, 1]. The operator denoted by c ∆ q is called the Caputo-type h-difference operator of order q for a function x : (hN) 0 → R and it is given by We recall the classical definition of the Z-transform associated to a sequence (see [25]). where z ∈ C is a complex number for which the series ∞ ∑ k=0 y(k)z −k converges absolutely.
The next proposition was proved in [26] and it will be useful in determining the characteristic equations that are associated to the considered system. Proposition 1. Let q ∈ (0, 1]. Defining y(j) = ( c ∆ q x)(t), where t ∈ (hN) 0 and t = jh, the Z-transform of y(j) is The following n-dimensional fractional-order system is considered: where q = (q 1 , q 2 , ..., q n ) ∈ (0, 1) n and f : (hN) 0 × R n → R n is a continuous function on the whole domain and Lipschitz-continuous with respect to the second variable, such that f (nh, 0) = 0 for any nh ∈ (hN) 0 . Let ϕ(nh, x 0 ) denote the unique solution of (1) satisfying the initial condition x(0) = x 0 . As in the case of the continuous counterpart explored in [27], we highlight the fact that, in general, because of the presence of the memory terms, the asymptotic stability of the null solution of system (1) will not be of exponential type [15,28]. Therefore, a particular type of non-exponential asymptotic stability, namely the notion of Mittag-Leffler stability was introduced for the case of fractional-order differential equations [29]. In the case of fractional-order difference systems, Mittag-Leffler stability was explored in [30]. In this paper, we focus our attention on a particular type of discrete Mittag-Leffler stability, specifically the O(n −q )-asymptotic stability, which pertains to the algebraic decay of the solutions. Therefore, we recall the following definitions of notions of stability: a. The trivial solution of (1) is called stable if for any ε > 0 there exists δ = δ(ε) > 0 such that for every x 0 ∈ R n satisfying x 0 < δ we have ϕ(nh, x 0 ) < ε for any n ≥ 0. b. The trivial solution of (1) is called asymptotically stable if it is stable and there exists ρ > 0, such that lim n→∞ ϕ(nh, x 0 ) = 0 whenever x 0 < ρ. c. The trivial solution of (1) is called O(n −q )-asymptotically stable if it is stable and there exists ρ > 0, such that for any x 0 < ρ one has ϕ(n, x 0 ) = O(n −q ) as n → ∞.

Stability Results for Systems of Two Fractional-Order Difference Equations
In what follows, we restrict our attention to the case of linear autonomous incommensurate fractional-order systems of the form: where the system's matrix A = (a ij ) is a real two-dimensional matrix, h is the discretization step, and q 1 , q 2 ∈ (0, 1) are the fractional orders of Caputo forward difference operators. By means of the Z-transform method, we obtain: where the Z-transforms of x and y are denoted by Z [x] = X and Z [y] = Y, respectively. This gives the following characteristic equation: which is equivalent to Observing that a 11 a 22 − a 12 a 21 = det A and multiplying the previous relation by h q 1 +q 2 , we obtain the characteristic function of the system (2): In the remainder of the paper, we will denote δ = det(A). Similarly to the continuous counterpart that was explored in [7,8], the next result can be obtained for the characterization of the stability and instability properties of system (2), with respect to the distribution of the roots of the characteristic function. 1. System (2) is O(n −q )-globally asymptotically stable if and only if all the roots of ∆ A (z) are inside the unit circle (|z| < 1), where q = min{q 1 , q 2 }. 2. If det(A) = 0 and ∆ A (z) has at least one root outside the closed unit circle (|z| ≥ 1), system (2) is unstable.

Fractional-Order Independent Results
We first explore several sufficient conditions that guarantee the instability of system (2), regardless of the choice of the fractional orders q 1 and q 2 .
Theorem 2 (Fractional-order independent instability results). System (2) is unstable for any choice of the fractional orders q 1 and q 2 if one of the following conditions hold: 3. a 11 > 0 and a 11 a 22 ≥ det(A) > 0.

Fractional-Order Dependent Results
In this section, let δ = det(A) > 0 and q 1 , q 2 ∈ (0, 1], h > 0 be arbitrarily fixed. In order to establish stability and instability conditions, the following curves and lines will be considered in the (a 11 , a 22 )-plane: We can easily prove the following lemmas, while using basic mathematical tools: Then ∆ A (−1) = 0 if and only if (a 11 , a 22 ) belongs to the line l(δ, q 1 , q 2 , h).

Proof. ∆ A has a root on the unit circle if and only if there exists
The following system is obtained by taking the real and imaginary parts of the terms in the above equation: (2 cos θ) q 1 +q 2 cos(q 1 + q 2 − 4)θ + a 11 h q 1 (2 cos θ) q 2 cos(q 2 − 2)θ + a 22 h q 2 (2 cos θ) q 1 cos(q 1 − 2)θ + h q 1 +q 2 δ = 0 (2 cos θ) q 1 +q 2 sin(q 1 + q 2 − 4)θ + a 11 h q 1 (2 cos θ) q 2 sin(q 2 − 2)θ + a 22 h q 2 (2 cos θ) q 1 sin(q 1 − 2)θ = 0 Solving this linear system with respect to the unknowns a 11 and a 22 we obtain: As in Lemma 2, ∆ A has a root on the unit circle of the form z = e i(π−2θ) = −e −2iθ , with θ ∈ 0, π 2 if and only if the following relation holds: By taking the real and imaginary parts of the previous equation, it follows that The second relation of system (5) is equivalent to which, combined with the first relation of (5) leads to Substituting θ in Equation (6), it follows that (a 11 , a 22 ) belongs to the line Λ(δ, q, h).

Lemma 4.
Let 0 < q 1 < q 2 < 1 and let a 11 (θ) and a 22 (θ) denote the functions that define the parametric equations of the curve Γ(δ, q 1 , q 2 , h). The following properties holds: Proof. The statements are proved in a straightforward manner.
Proof of Statement 1. Using the L'Hospital rule, we obtain: Hence: It follows that (a 0 11 , a 0 22 ) ∈ l(δ, q 1 , q 2 , h). Proof of Statement 2. The proof is trivial and will be omitted; it follows from the fact that lim θ→ π Proof of Statement 3. Dropping for simplicity the arguments of the functions ρ i (q 1 , q 2 , θ) and u(θ, h) and applying the inequality of arithmetic and geometric means, we obtain: as one can show, by elementary trigonometric inequalities, that We have proved in Lemma 2 that, if (a 11 , a 22 ) ∈ Γ(δ, q 1 , q 2 , h), then the characteristic function ∆ A admits a pair of complex conjugated roots on the unit circle and, therefore, in line with the classical theory of discrete dynamical systems, a discrete Hopf (Neimark-Sacker) bifurcation is expected to take place. Moreover, if (a 11 , a 22 ) ∈ l(δ, q 1 , q 2 , h), by Lemma 1, the characteristic function ∆ A has a root z = −1, which is normally associated with a Flip bifurcation in the case of discrete dynamical systems (see [31]). Still, it is of high importance to underline that, for systems of fractional-order difference equations, there are no thorough results regarding their bifurcation theory, which remains yet to be investigated.

Lemma 5.
Let α ∈ R and q ∈ (0, 1). All of the roots of the equation are in the interior of the unit disc if and only if α ∈ (0, 2 q ).
By the implicit function theorem, taking into consideration the characteristic Equation (3), we deduce: where P(z) = ∂ ∂z ∆ A (z). On one hand, we have: On the other hand, it can be shown that We obtain: Therefore, this leads to: A laborious but straightforward calculation (which can also by verified by symbolic computation tools) gives ∂|z| In a similar manner, we compute ∂|z| ∂a 22 Taking into account the parametric equations of the curve Γ(δ, q 1 , q 2 , h), it follows that the gradient vector ∇|z|(a * 11 , a * 22 ) is normal to the curve Γ(γ, q 1 , q 2 , q), pointing towards the region above the curve. Hence, the following transversality condition can be expressed in terms of the directional derivative: ∇ u |z|(a * 11 , a * 22 ) = ∇|z|(a * 11 , a * 22 ), u > 0, for any vector u pointing towards the region above the curve Γ(δ, q 1 , q 2 , h). Therefore, as the parameters (a 11 , a 22 ) vary and cross the curve Γ(δ, q 1 , q 2 , h) into the region above the curve, |z| becomes larger than 1. Case 2. If (a * 11 , a * 22 ) ∈ l(δ; q 1 , q 2 , h), based on Lemma 1, we will now denote, by z(a 11 , a 22 , δ, q 1 , q 2 ), the root of the characteristic function ∆ A satisfying the relation z(a * 11 , a * 22 , δ, q 1 , q 2 ) = −1.
Following the same method as above, Equation (10) now provides: As (a * 11 , a * 22 ) ∈ l(δ; q 1 , q 2 , h), it follows that Hence: In a similar manner, computing ∂|z| ∂a 22 It is easy to see that the vector (h/2) q 1 , (h/2) q 2 is normal to the line l(δ, q 1 , q 2 , h). Moreover, as a * 11 < a 0 11 , i.e., the gradient vector ∇|z|(a * 11 , a * 22 ) is normal to the line l(γ, q 1 , q 2 , h), pointing towards the region that is below the line. Similarly as in the previous case, the transversality condition implies that, as the parameters (a 11 , a 22 ) vary and cross the line l(δ, q 1 , q 2 , h) into the region below the line, |z| becomes larger than 1.
Therefore, from the previously discussed case, we deduce that, as the parameters (a 11 , a 22 ) are varied and leave the domain S(δ, q 1 , q 2 , h), the asymptotic stability of the system (2) is lost.
Proof of statement 2. The proof follows the same lines as previously and, hence, it will be omitted.
We remark that inequality (8) plays an important role in the proof of Theorem 3, giving an upper bound for the discretization step h. Moreover, Figure 1 exemplifies the results that were obtained in Theorem 3, for several values of the discretization step h, when considering the fractional orders q 1 = 0.8 and q 2 = 0.4 and δ = det(A) = 5.

A Discrete FitzHugh-Nagumo Neuronal Model
The FitzHugh-Nagumo neuronal model describes a biological neuron's spiking behavior and it is a simplification of the well-known Hodgkin-Huxley model. Several studies have previously suggested that, at the neuronal membrane level, conductance adaptation is history dependent and, hence, its response can be described by a power-law function [21,22,32], leading to mathematical models that involve fractional-order differential equations [33]. Here, we consider a discretization of the fractional-order version of the FitzHugh-Nagumo model that was studied in [8]: where v is the neuronal membrane potential, w represents a recovery variable, and I is an external excitation current. Moreover, q 1 and q 2 are the fractional orders of the Caputo h-difference operators, with 0 ≤ q 1 ≤ q 2 ≤ 1.
It is important to emphasize that, in the numerical investigation of the continuous-time fractional-order FitzHugh-Nagumo model undertaken in [8], the well-known Adams-Bashforth-Moulton scheme was used [34]. Instead, in this section, we investigate the discrete-time model that was obtained by replacing the fractional-order Caputo derivatives with fractional-order h-differences of Caputo type.
Rewriting the second equation of system (11), we obtain: where α = 1 d > 1, β = c d and φ = rd ∈ (0, 1). It follows that system (11) is a conductance-based system of the form: with and w ∞ (v) = αv + β is a linear function. Similarly to the case of the continous version of the fractional-order FitzHugh-Nagumo neuronal model [8], the equilibrium states of the fractional-order neuronal model (12) can be determined. In order to investigate the stability of equilibrium states, the Jacobian matrix associated to system (12) at an arbitrary equilibrium state (v * , w * ) = (v * , w ∞ (v * )) is: Based on the results that were obtained in Theorem 3, the stability regions have been determined and transposed to the (q 1 , q 2 )-plane for different values of the equilibrium membrane potential v * and several discretization steps. Figure 2 plots these regions.  Taking r = 0.08, c = 0.7, d = 0.8, and I = 1.25 as values for the parameters of system (11), the equilibrium potential is v * = 0.804848. In Figure 3, the evolution of the membrane potential can be observed, taking an initial condition from a small neighborhood of the equilibrium. We considered the value of the fractional order for the recovery variable to be q 2 = 1 and varied the value of the fractional order corresponding to the membrane potential q 1 . A Hopf bifurcation appears to be occurring for a value q 1 in the interval (0.47, 0.48) as for q 1 = 0.47, we have an asymptotically stable behaviour and for q 1 ≥ 0.48, oscillations can be observed, depicting individual spikes that can be modulated by the fractional order and that complement the dynamic behavior displayed by classical integer order counterpart, as observable in the last image of Figure 3.

Conclusions
Theoretical results regarding both asymptotic stability and instability properties of two-dimensional systems of autonomous linear fractional-order difference equations of Caputo type have been explored. The theoretical findings were later applied in the investigation of a fractional-order version of the well-known FitzHugh-Nagomo neuronal model. Moreover, numerical simulations revealed that the chaotic bursting and spiking behavior of the membrane potential of a biological neuron can be modulated by varying the corresponding fractional order of the system. This is in accordance with previous results that were published in the literature for other types of fractional-order neuronal models [22,[35][36][37][38], where the spike frequency and amplitude reduction have also been observed as the fractional-order decreases (which is associated to a greater influence of the memory of the membrane potential).