**2. Fractional-Order Glucose-Insulin Regulatory System**

Because diverse biological systems present memory effects, direct generalization maybe by applying fractional-order differential equations, i.e., arbitrary (non-integer) orders. In this manner, we can get a closer insight into the real phenomena. The benefits of fractional-order are mainly to capture the entire time evolution for physical processes, being a kind of memory index, as well as more degrees of freedom for the resulting models. The analysis conducted in this work is inspired by the glucose-insulin regulatory system in [46], as given in Equation (1). This model was derived from the Ackerman, Bajaj, and Rao, and predator-prey Volterra models. As a consequence, the resulting system is formulated as the trade-off relating the glucose and insulin, but including the beta-cells interaction. Apart from normal metabolic conditions described by the classical relation between prey (glucose) and predator (insulin), the model (1) also characterizes the metabolic disorders as chaotic oscillations. Abnormal biological conditions, including glucose-insulin interactions, can be stated as chaotic evolutions of the dynamical systems, as given in [47–50], to mention a few. From experimental data, Ref. [51] also discussed that the T1DM may have chaotic behaviors. In this manner, the proposed analysis of the glucose-insulin regulatory system is in the same line of previous studies, but instead applying fractional-calculus theory.

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

we now introduce a fractional-order version [52], but instead applying the Caputo's definition. Let us consider the fractional differential operator *<sup>D</sup><sup>q</sup> <sup>f</sup>*(*t*) = *<sup>J</sup>n*−*γD<sup>n</sup> <sup>f</sup>*(*t*), with *<sup>q</sup>* <sup>&</sup>gt; 0 and *<sup>n</sup>* <sup>∈</sup> <sup>N</sup>, where *D<sup>n</sup>* describes the *n*-order derivative, while *J<sup>γ</sup>* defines the *γ*-order integral operator in the sense of Riemann-Liouville as

**Definition 1.** *For a function f*(*t*)*, the fractional-order (γ) integral, with <sup>γ</sup>* <sup>∈</sup> <sup>R</sup><sup>+</sup> *can be determined by*

$$J^{\gamma}f(t) = \frac{1}{\Gamma(\gamma)} \int\_{t\_0}^{t} (t-s)^{\gamma-1} f(s) ds,\tag{2}$$

*where t* <sup>≥</sup> *<sup>t</sup>*<sup>0</sup> *and* <sup>Γ</sup>(·) *is the Gamma function, defined as* <sup>Γ</sup>(*s*) = <sup>∞</sup> <sup>0</sup> *t <sup>s</sup>*−1*e*−*<sup>t</sup> dt.*

In this manner, the Caputo operator of fractional-order can be established by

**Definition 2.** *The fractional derivative of Caputo with order <sup>q</sup> for a function <sup>f</sup>*(*t*) <sup>∈</sup> <sup>C</sup>*n*([*t*0, <sup>∞</sup>), <sup>R</sup>) *is given as*

$$D^q f(t) = \frac{1}{\Gamma(n-q)} \int\_{t\_0}^t \frac{f^{(n)}(s)}{(t-s)^{q-n+1}} ds,\tag{3}$$

*and when* 0 < *q* < 1 *is given, as follows*

$$D^q f(t) = \frac{1}{\Gamma(1-q)} \int\_{t\_0}^t \frac{f'(s)}{(t-s)^q} ds. \tag{4}$$

By substituting Equation (4) in Equation (1), we get the fractional-order glucose-insulin regulatory system that is given by

$$\begin{array}{rcl}D^{\theta\_1}\mathbf{x} &=& -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\tau\mathbf{z}^3 - 0.19, \\ D^{\theta\_2}\mathbf{y} &=& -a\_8\mathbf{x}\mathbf{y} + 3.84\mathbf{x}^2 + 1.2\mathbf{x}^3 + 0.3y(1 - \mathbf{y}) - 1.37\mathbf{z} + 0.3\mathbf{z}^2 - 0.22\mathbf{z}^3 - 0.56, \\ D^{\theta\_3}\mathbf{z} &=& 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{5}$$

where *Dqi*(*x*, *y*, *z*) with *i* = 1, 2, 3 is the the fractional derivative operator in Caputo's sense [53], *qi* is the fractional-order satisfying 0 < *qi* ≤ 1, (*x*, *y*, *z*) describe the population density of insulin, the population size of glucose, and the population size of *β*-cells, respectively. The representative reduction of insulin levels for glucose deficiency is described by *a*1, whereas the parameter *a*<sup>8</sup> stands for the impact of insulin on glucose. Additionally, the increment rate of insulin concentrations released by b-cells is provided by *a*7. Finally, the increasing rate of the *β*-cells provoked, when the glucose levels also grow, is given by *a*<sup>15</sup> [46]. From the nature of how this system was conceived, the parameters must be nonnegative [46].

The discrete-time version for the proposed system (5) is set as:

$$\begin{aligned} x\_{n+1} &= x\_0 + \frac{h^{q\_1}}{\Gamma(q\_1 + 2)} \left( f\_1(x\_{n+1}^p, y\_{n+1}^p, z\_{n+1}^p) + \sum\_{k=0}^n a\_{1,k,n+1} f\_1(x\_k, y\_k, z\_k) \right), \\ y\_{n+1} &= y\_0 + \frac{h^{q\_2}}{\Gamma(q\_2 + 2)} \left( f\_2(x\_{n+1}^p, y\_{n+1}^p, z\_{n+1}^p) + \sum\_{k=0}^n a\_{2,k,n+1} f\_2(x\_k, y\_k, z\_k) \right), \\ z\_{n+1} &= z\_0 + \frac{h^{q\_3}}{\Gamma(q\_3 + 2)} \left( f\_3(x\_{n+1}^p, y\_{n+1}^p, z\_{n+1}^p) + \sum\_{k=0}^n a\_{3,k,n+1} f\_3(x\_k, y\_k, z\_k) \right), \end{aligned} \tag{6}$$

where

$$\begin{aligned} x\_{n+1}^p &= x\_0 + \frac{1}{\Gamma(q\_1)} \sum\_{k=0}^n \beta\_{1,k,n+1} f\_1(\mathbf{x}\_k, y\_k, z\_k), \\ y\_{n+1}^p &= y\_0 + \frac{1}{\Gamma(q\_2)} \sum\_{k=0}^n \beta\_{2,k,n+1} f\_2(\mathbf{x}\_k, y\_k, z\_k), \\ z\_{n+1}^p &= z\_0 + \frac{1}{\Gamma(q\_3)} \sum\_{k=0}^n \beta\_{3,k,n+1} f\_3(\mathbf{x}\_k, y\_k, z\_k), \end{aligned} \tag{7}$$

with

$$a\_{i,k,n+1} = \begin{cases} n^{q\_i+1} - (n - q\_i)(n+1)^{q\_i}, k = 0, & k = 0\\ (n - k + 2)^{q\_i+1} + (n - k)^{q\_i+1} - 2(n - k + 1)^{q+1}, & 1 \le k \le n, \\ 1, & k = n+1, \end{cases} \tag{8}$$

and

$$
\beta\_{i,k,n+1} \frac{h^{q\_i}}{q\_i} \left( (n+1-k)\_i^q - (n-k)\_i^q \right). \tag{9}
$$

note that *<sup>x</sup>*0, *<sup>y</sup>*0, *<sup>z</sup>*<sup>0</sup> are the initial values for, *<sup>f</sup>*1(*x*, *<sup>y</sup>*, *<sup>z</sup>*) = −*a*1*<sup>x</sup>* + 0.1*xy* + 1.09*y*<sup>2</sup> − 1.08*y*<sup>3</sup> + 0.03*<sup>z</sup>* − 0.06*z*<sup>2</sup> + *<sup>a</sup>*7*z*<sup>3</sup> − 0.19, *<sup>f</sup>*2(*x*, *<sup>y</sup>*, *<sup>z</sup>*) = −*a*8*xy* + 3.84*x*<sup>2</sup> + 1.2*x*<sup>3</sup> + 0.3*y*(<sup>1</sup> − *<sup>y</sup>*) − 1.37*<sup>z</sup>* + 0.3*z*<sup>2</sup> − 0.22*z*<sup>3</sup> − 0.56 and *<sup>f</sup>*3(*x*, *<sup>y</sup>*, *<sup>z</sup>*) = *<sup>a</sup>*15*<sup>y</sup>* − 1.35*y*<sup>2</sup> + 0.5*y*<sup>3</sup> + 0.42*<sup>z</sup>* + 0.15*yz*, whether *qi* = *<sup>q</sup>*, *<sup>i</sup>* = 1, 2, 3 the system is said to be commensurate then the convergence order is described as <sup>|</sup>*y*(*tn*) <sup>−</sup> *yn*<sup>|</sup> <sup>=</sup> *<sup>O</sup>*(*h*min(2,1+*q*)), *h* → 0 [54–56].
