**3. Stability Analysis of Fractional-Order Glucose-Insulin Model**

A non-local differential equation with fractional-order *q* ∈ (0, 1), normally has a stability region larger than that of the integer-order version with *q* = 1 [57]. We introduce the following theorems and definitions to discuss the stability of the proposed fractional-order system.

**Definition 3.** *A system denoted as*

$$D^{\eta\_i} \mathbf{x}\_i(t) = f\_i(\mathbf{x}\_1(t), \mathbf{x}\_2(t), \dots, \mathbf{x}\_n(t), t), \quad \mathbf{x}\_i(0) = \mathbf{c}\_i, \quad i = 1, 2, \dots, n,\tag{10}$$

*with trajectory x*(*t*) = 0*, where ci are the starting conditions is t* <sup>−</sup>*<sup>q</sup> asymptotically stable if there is a nonnegative real q, so that:*

> ∀||*x*(*t*)||*with t* <sup>≤</sup> *<sup>t</sup>*0, <sup>∃</sup>*N*(*x*(*t*)),*such that* <sup>∀</sup>*<sup>t</sup>* <sup>&</sup>gt; *<sup>t</sup>*0, ||*x*(*t*)|| ≤ *Nt*−*<sup>q</sup>* .

*Subsequently, fractional-order systems are known to have long memory, since they obey a behavior like t* <sup>−</sup>*q, which slowly tends to 0 for the solutions x(t). A particular case of Mittag–Leffler stability is stated as power-law stability t*<sup>−</sup>*q. [58,59].*

**Theorem 1.** *Let a commensurate order system described by <sup>D</sup>qx* <sup>=</sup> *Ax, with <sup>x</sup>*(0) = *<sup>x</sup>*0*,* <sup>0</sup> <sup>&</sup>lt; *<sup>q</sup>* <sup>&</sup>lt; <sup>1</sup>*, <sup>x</sup>* <sup>∈</sup> <sup>R</sup>*<sup>n</sup> and <sup>A</sup>* <sup>∈</sup> <sup>R</sup>*n*×*n, is asymptotically stable if* <sup>|</sup> arg(*λ*)<sup>|</sup> <sup>&</sup>gt; *<sup>q</sup>π*/2, *is fulfilled for all eigenvalues of A. Besides, the critical eigenvalues fulfilling* | arg(*λ*)| = *qπ*/2 *holding geometric multiplicity of one, where the geometric multiplicity of an eigenvalue is called the dimension of the subspace vectors v for Av* = *λv [60–65].*

**Theorem 2.** *Let <sup>D</sup>qx* = *Ax*, *with <sup>x</sup>*(0) = *<sup>x</sup>*<sup>0</sup> *an incommensurate-order system, where <sup>x</sup>* = (*x*1, *<sup>x</sup>*2, ... , *xn*)*<sup>T</sup>* ∈ <sup>R</sup>*n, <sup>D</sup>qx* = (*Dq*<sup>1</sup> *<sup>x</sup>*1, *<sup>D</sup>q*<sup>2</sup> *<sup>x</sup>*2, ... , *<sup>D</sup>qn xn*)*T*, *qi* <sup>∈</sup> <sup>R</sup>+*, <sup>i</sup>* <sup>=</sup> 1, 2, . . . , *<sup>n</sup>*.*, qi* <sup>∈</sup> (0, 1)*, and <sup>A</sup>* = (*aij*) <sup>∈</sup> <sup>R</sup>*n*×*n, i* = *j* = 1, 2, . . . , *n*. *By supposing w as the lowest common multiple of the denominators ui's of qi's, with qi* = *vi*/*ui,* (*ui*, *vi*) = <sup>1</sup>*, ui*, *vi* <sup>∈</sup> <sup>Z</sup><sup>+</sup> *for i* <sup>=</sup> 1, 2, . . . , *<sup>n</sup>*.*, the matrix of system is given by*

$$
\Delta(\lambda) = \begin{bmatrix}
\lambda^{\text{u}\eta\_1} - a\_{11} & -a\_{12} & \dots & -a\_{1n} \\
\vdots & \vdots & \ddots & \vdots \\
\end{bmatrix} \tag{11}
$$

*Therefore, the system Dqx* = *Ax is globally asymptotically stable if the roots λ's of its characteristic equation* det(Δ(*λ*)) = 0 *fulfil* | arg(*λ*)| > *π*/2*w, [60–65].*

**Theorem 3.** *The equilibria E*<sup>∗</sup> *with a instability measure defined as*

$$\rho = (\pi/2w) - \min\_{i} \{ |\arg(\lambda\_i)| \},\tag{12}$$

*is asymptotically stable if and only if* (12) *is defined nonpositive. Where λi's are the roots of characteristic equation:* det(*diag*([*λwq*<sup>1</sup> *<sup>λ</sup>wq*<sup>2</sup> ... *<sup>λ</sup>wqn* ]) <sup>−</sup> *<sup>∂</sup> <sup>f</sup>* /*∂x*|*x*=*E*<sup>∗</sup> ) = <sup>0</sup>*,* <sup>∀</sup>*E*<sup>∗</sup> <sup>∈</sup> <sup>Ω</sup> *[64,65].*

**Remark 1.** *If ρ is nonnegative, the underlying system can show a chaotic evolution since the equilibrium E*<sup>∗</sup> *is unstable [64,65].*

*Symmetry* **2020**, *12*, 1395

The first step consists of finding the equilibrium points of (5), which are calculated via **f**(**x**) = 0, as follows

$$\begin{array}{llll} 0 & = & -a\_1 \mathbf{x} + 0.1 \mathbf{x} \mathbf{y} + 1.09 \mathbf{y}^2 - 1.08 \mathbf{y}^3 + 0.03 \mathbf{z} - 0.06 \mathbf{z}^2 + a\_7 \mathbf{z}^3 - 0.19, \\ 0 & = & -a\_8 \mathbf{x} \mathbf{y} + 3.84 \mathbf{x}^2 + 1.2 \mathbf{z}^3 + 0.3 \mathbf{y} (1 - \mathbf{y}) - 1.37 \mathbf{z} + 0.3 \mathbf{z}^2 - 0.22 \mathbf{z}^3 - 0.56, \\ 0 & = & a\_{15} \mathbf{y} - 1.35 \mathbf{y}^2 + 0.5 \mathbf{y}^3 + 0.42 \mathbf{z} + 0.15 \mathbf{y} \mathbf{z}. \end{array} \tag{13}$$

Because of the biological interpretation of state variables [46], the stability analysis for *E*∗ = (*x*∗, *y*∗, *z*∗) is only performed for nonnegative fixed points. Equation (5) has the following Jacobian

$$J = \begin{pmatrix} -a\_1 + 0.1y\* & 0.1x\* + 2.18y^\* - 3.24y^{\*2} & 0.03 - 0.12z^\* + 3a\gamma z^{\*2} \\ -a\_8y^\* + 7.68x^\* + 3.6x^{\*2} & -a\_8x + 0.3(1 - 2y^\*) & -1.37 + 0.6z - 0.66z^{\*2} \\ 0 & a\_{15} - 2.7y^\* + 1.5y^{\*2} + 0.15z & 0.42 + 0.15y^\* \end{pmatrix}. \tag{14}$$

By setting *a*<sup>1</sup> = 1.3, *a*<sup>7</sup> = 2.01, *a*<sup>8</sup> = 0.22, *a*<sup>15</sup> = 0.3, we compute the equilibrium points and eigenvalues, as shown in Table 1. As can be noted, the system (5) has two nonnegative fixed points characterized for saddle points of index-1 and index-2, respectively. According to Theorem 1, chaos behavior may arise in the fractional-order glucose-insulin system when | arg(*λ*)| > *qπ*/2, holds. Therefore, the minimum commensurate fractional-order leading to chaotic oscillations is *q* > 0.7638, which is remarkably lower than that reported in [52].

**Table 1.** Nonnegative equilibrium points and eigenvalues of the system (5).


The chaotic attractor existence and local stability of the equilibria are presented below. The Jacobian matrix considering the equilibria *E*<sup>1</sup> and parameter *a*<sup>1</sup> is

$$J\_{E\_1} = \begin{pmatrix} -a\_1 + 0.1853 & -7.0110 & 9.8541 \\ 8.0788 & -0.9887 & -1.6903 \\ 0 & 0.6420 & 0.6420 \end{pmatrix},\tag{15}$$

which generates the characteristic equation given, as follows

$$P(\lambda) = \lambda^3 + (a\_1 + b\_1)\lambda^2 + (b\_2a\_1 + b\_3)\lambda + b\_4a\_1 - b\_5 = 0,\tag{16}$$

where *b*<sup>1</sup> = 0.1054, *b*<sup>2</sup> = 0.2907, *b*<sup>3</sup> = 56.9826, *b*<sup>4</sup> = 0.395, and *b*<sup>5</sup> = 90.7232, when considering the Routh–Hurwitz theory, obtaining the following conditions *a*<sup>1</sup> + *b*<sup>1</sup> > 0, *b*2*a*<sup>1</sup> + *b*<sup>3</sup> > 0, and *b*4*a*<sup>1</sup> − *b*<sup>5</sup> > 0, which become in *a*<sup>1</sup> > −*b*<sup>1</sup> and *b*<sup>3</sup> > −*b*2*a*1, respectively. Due *b*<sup>5</sup> > *b*4*a*1, there is a root positive which is unstable and a pair of complex conjugate with the negative real part. Thus *E*<sup>1</sup> is unstable being a saddle point of index 1.

The presence of a Hopf bifurcation in the system (5) at the equilibrium point *E*<sup>1</sup> is analyzed by replacing *λ* = *iω* in (16) obtaining

$$(i\dot{w})^3 + (a\_1 + b\_1)(i\omega)^2 + (b\_2a\_1 + b\_3)(i\omega) + b\_4a\_1 - b\_5 = 0,\tag{17}$$

then

$$-\left(iw\right)^{3} - \left(a\_{1} + b\_{1}\right)\omega^{2} + \left(b\_{2}a\_{1} + b\_{3}\right)\left(i\omega\right) + b\_{4}a\_{1} - b\_{5} = 0,\tag{18}$$

taking into account the real part of (18), therefore,

$$
\omega^2 = \frac{b\_5 - b\_4 a\_1}{-(a\_1 + b\_1)},
\tag{19}
$$

because *a*<sup>1</sup> + *b*<sup>1</sup> > 0, and *b*5, *b*4, *a*<sup>1</sup> > 0, hence *ω*<sup>2</sup> < 0, so it is not possible, then, the equilibrium point *E*<sup>1</sup> does not suffer a Hopf bifurcation.

In order to study the stability of *E*<sup>1</sup> we consider the discriminant *D*(*P*) of characteristic equation *P*(*λ*), given as follows

$$D(P) = 18d\_1d\_2d\_3 + (d\_1d\_2)^2 - 4d\_3d\_1^3 - 4d\_2^3 - 27d\_3^2. \tag{20}$$

where *d*<sup>1</sup> = (*a*<sup>1</sup> + *b*1), *d*<sup>2</sup> = (*b*2*a*<sup>1</sup> + *b*3), and *d*<sup>3</sup> = (*b*4*a*<sup>1</sup> − *b*5). Following the Routh–Hurwitz stability conditions for fractional-order differential equations [66], we establish the following terms:


**Remark 2.** *For parameter values di, i* = 1, 2, 3 *with d*<sup>1</sup> = 1.4054*, d*<sup>2</sup> = 57.3606 *and d*<sup>3</sup> = −90.2097 *the discriminant <sup>D</sup>*(*P*) = −1.098041820 × <sup>10</sup>6*, but <sup>d</sup>*1*d*<sup>2</sup> = 80.6144 = *<sup>d</sup>*3*, and <sup>d</sup>*<sup>3</sup> < <sup>0</sup> *thence, the Routh–Hurwitz conditions are unsatisfied. Thus, the E*<sup>1</sup> *is unstable for the given parameters.*

When *q*<sup>1</sup> = *q*<sup>2</sup> = *q*<sup>3</sup> ≡ *q* = 0.9 = 9/10, with *w* = 10, the characteristic equation of (5) in the equilibrium point *E*<sup>1</sup> = (0.802, 1.853, 1.286) is given by Theorem 3, as follows

$$
\lambda^{27} + (a\_1 + b\_1)\lambda^{18} + (b\_2a\_1 + b\_3)\lambda^9 + b\_4a\_1 - b\_5 = 0,\tag{21}
$$

if we set *a*<sup>1</sup> = 1.3 the characteristic equation has an unstable root *λ* = 1.0434 and the | arg(*λ*)| < *π*/20, moreover the *E*<sup>1</sup> is saddle point. While the characteristic equation at the equilibrium point *E*<sup>2</sup> = (0.584, 0.832, 0.728) is

$$
\lambda^{27} + (a\_1 - 0.3002)\lambda^{18} + (-0.217a\_1 + 0.868)\lambda^9 - 1.2036a\_1 + 12.88179578 = 0,\tag{22}
$$

by considering *a*<sup>1</sup> = 1.3 we obtain unstable roots *λ*1,2 = 1.0771±0.1444*i*. thence, the instability measure of the system is *ρ* = (*π*/2*w*) − 0.1333 > 0. Therefore, the fractional-order system (5) satisfies the necessary condition for exhibiting a chaotic attractor according to Theorem 3. This can be understood by locating the respective eigenvalues in the complex plane. Figure 1, displays the 27 eigenvalues of the system in the complex plane, the unstable region is delimited by the red lines which is denoted by *π*/2*w*. The eigenvalues for *E*<sup>1</sup> are given in Figure 1a, we can observe that the system has one eigenvalue in the unstable region, it is a saddle point of index 1. Meanwhile, Figure 1b shows the eigenvalues of equilibrium point *E*2; it has two unstable eigenvalues (saddle point of index 2).

**Figure 1.** Eigenvalues of the system (5) in the complex plane. (**a**) eigenvalues for equilibria *E*1, (**b**) eigenvalues for equilibria *E*2.

Figure 2(a)-(c), exhibits the chaotic attractor found in the glucose-insulin model (5) for fractional-orders *q*<sup>1</sup> = *q*<sup>2</sup> = *q*<sup>3</sup> = 0.9. The results were computed by applying the predictor-corrector scheme Adams–Bashforth–Moulton (ABM) [54–56]. Due to Caputo's fractional differential operator (4) permits to select both homogeneous and inhomogeneous initial conditions, the ABM algorithm can be executed without particular constraints [56].

As well known, the differential equations that model dynamical systems regularly present symmetries and, therefore, it is rare the solutions do not evolve in a symmetrical orientation. The most common effect is asymmetric attractors for a proper range of parameters, which almost all converges to a singlesymmetric attractor [67]. For dynamical systems with three state variables *x*, *y*, and *z*, there are three types of involuntary symmetries: inversion, rotation, and reflection. We consider the following transformations (*x*, *y*, *z*) → (−*x*, −*y*, −*z*), (*x*, *y*, *z*) → (−*x*, −*y*, *z*), and (*x*, *y*, *z*) → (−*x*, *y*, *z*) corresponding to the invariance of the equations with changes of sign in three, two, and one variable, respectively. As a result, the system (5) does not have a symmetry under the proposed transformations, and its trajectories (*x*(*t*), *y*(*t*), *z*(*t*)) cannot cross the (0, 0, 0) coordinate.

**Figure 2.** Chaotic attractor of the fractional-order glucose-insulin model (5) for fractional-orders *q* = 0.9, *a*<sup>1</sup> = 1.3, *a*<sup>7</sup> = 2.01, *a*<sup>8</sup> = 0.22, *a*<sup>15</sup> = 0.3, (*x*0, *y*0, *z*0)=(0.5, 1.2, 1), and integration step-size *h* = 0.01. (**a**) *x* − *y* phase plane, (**b**) *x* − *z* phase plane, (**c**) *y* − *z* phase plane.
