EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 2, Article Number 6152 ISSN 1307-5543 – ejpam.com Published by New York Business Global Hyers-Ulam Stability and Control of Fractional Glucose-Insulin Systems Sayed Saber1,2,∗, Brahim Dridi3, Abdullah Alahmari3, Mohammed Messaoudi4 1 Mathematics Department, Faculty of Science, Al-Baha University, Saudi Arabia 2 Department of Mathematics and Computer Science, Faculty of Science, Beni-Suef University, Egypt 3 Mathematics Department, Faculty of Sciences, Umm Al-Qura University, P.O. Box 14035, Makkah, 21955, Saudi Arabia 4 Imam Mohammad Ibn Saud Islamic University (IMSIU), College of Science, Department of Mathematics and Statistics, Riyadh, Saudi Arabia Abstract. This paper presents a novel fractional-order model for glucose-insulin dynamics using the Caputo-Fabrizio (CF) derivative, which accounts for memory effects through its nonsingular exponential kernel. Existence and uniqueness of the solution are established via fixed point theory, and infinite series solutions are obtained using the Sumudu transform. Hyers-Ulam stability is analyzed to assess the system’s robustness against perturbations. A linear control strategy is introduced to regulate glucose levels, demonstrating potential for integration with real-time insulin delivery systems. Compared to classical integer-order models, the proposed approach provides improved accuracy, enhanced stability, and deeper insight into the chaotic behavior of glucose- insulin interactions. This framework supports the development of personalized diabetes treatment and adaptive control strategies. 2020 Mathematics Subject Classifications: 34L99 Key Words and Phrases: Fractional calculus, Glucose-insulin model, Chaos control, Stability analysis, Caputo fractional derivative, Numerical simulation, Diabetes modeling. 1. Introduction The prevalence of diabetes among adults worldwide has increased to over 422 million, according to WHO [1]. Diabetes and its complications have been studied extensively over the last decade, highlighting the need for long-term management. Diabetes burden can be reduced through early detection and prevention [2]. A healthcare professional is essential ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.6152 Email addresses: Sayed011258@science.bsu.edu.eg (S. Saber), iodridi@uqu.edu.sa (B. Dridi), aaahmari@uqu.edu.sa (A. Alahmari), mmessaoudi@imamu.edu.sa (M. Messaoudi) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 2 of 24 in managing a patient’s condition and preventing complications [3]. Literature on glucose- insulin dynamics, epidemiology, complications, and the economic burden of diabetes is growing. Many mathematical, statistical, and computational algorithms have been pro- posed to help us understand diabetes [4, 5]. The studies cover the range of topics discussed above, such as the dynamics of glucose and insulin in relation to insulin [6, 7], computer algorithms and devices [8–10], and the costs associated with diabetes [11, 12]. Addition- ally, research has been conducted on the psychosocial aspects of diabetes, such as quality of life and social support [13, 14], and its impact on public health [15]. A recent develop- ment in fractional calculus has been its ability to model complex systems with arbitrary differentiation orders. Diabetes research has benefited from enhanced modeling of glucose- insulin dynamics, as well as complications associated with diabetes [16, 17]. In addition to the Riemann-Liouville, Caputo, and Caputo-Fabrizio (CF) fractional calculus operators, Atangana-Baleanu can be used to capture power laws and exponential decay [18]. The CF operator, which has a nonsingular kernel, has successfully been applied to glucose-insulin modeling, providing insights into memory effects and chaos within the system [19–54]. Sadek et al. (2024) introduced novel Θ-fractional operators, expanding the framework of fractional calculus and demonstrating their efficiency in optimization problems using the Galerkin-Bell method [55]. In [56], Sadek proposed a cotangent fractional derivative, providing an alternative approach to solving complex differential equations with applica- tions in fractional modeling. In the context of epidemiology, Sadek et al. [57] developed a fractional-order model to predict COVID-19 dynamics in Morocco, incorporating isola- tion and vaccination strategies, where the inclusion of memory effects improved predictive accuracy [58]. Additionally, Sadek et al. (2022) applied fractional calculus to model the synthesis of TiO2 nanopowder via the sol–gel method at low temperatures, capturing the complexities of nanopowder formation and improving process modeling in nanotechnol- ogy [59]. These contributions collectively enhance the theoretical and applied aspects of fractional calculus, demonstrating its versatility across diverse scientific domains. An exponential decay function is used to model the fractional differential operator to understand how memory impacts fractional glucose-insulin regulation. With this ker- nel, we can study memory effects. Using Caputo-Fabrizio fractional operators to analyze fractional variable-order glucose-insulin regulation, memory effects enhance stability and synchronization. It highlights the potential for developing more effective diabetes con- trol strategies based on uniqueness and boundedness. We also examine Caputo-Fabrizio fractional variable-order glucose-insulin regulation, demonstrating that fractional chaotic systems can be synchronized. Fixed point theory is used to prove uniqueness and bound- edness of solutions, and CF fractions are used to analyze Hyers-Ulam stability. Sumudu transforms are used to derive infinite series solutions that demonstrate efficiency. Different smooth functions defined in (0,1] can be used to simulate fractional derivatives. Our novel fractional-order glucose-insulin model incorporates memory effects and non- local interactions based on the Caputo-Fabrizio derivative. Compared to traditional integer-order models, this approach offers a more accurate representation of physiolog- ical processes. Using Sumudu transforms for solution derivation enhances computational efficiency and reduces numerical errors. Based on Hyers-Ulam stability analysis, our model S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 3 of 24 is robust against perturbations, while a newly developed linear control strategy effectively regulates glucose levels, with potential applications in real-time insulin delivery. Based on the results, fractional-order derivatives provide a more accurate and stable representation of glucose-insulin fluctuations. Moreover, our model simulates disease progression, opti- mizes insulin dosage for personalized treatment, and stabilizes extreme glucose variations, which prevents severe complications such as hypoglycemia and hyperglycemia. In contrast to previous studies, our approach provides novel insights into chaotic dynamics in glucose regulation, which contributes to the efficient management of diabetes. A key aspect of this paper is the utilization of Caputo-Fabrizio-Caputo fractional derivatives, which possess a non-singular kernel, enabling an accurate description of di- verse physiological processes. Diabetes patients can analyze glucose-insulin interactions and gain insight into memory and chaotic behaviors with the CF-based glucose-insulin regulation model. The model can also be used to develop advanced control strategies. The proposed model outperforms glucose-insulin models with integer-order derivatives. Real-world implementation of this model may be challenging due to fractional deriva- tives and variable orders. Additionally, collecting and analyzing physiological parameters is complex. The model must be validated across diverse patient populations to ensure robustness and reliability. The CF fractional-order glucose-insulin system requires linear controllers to achieve equilibrium. Diabetes management is significantly improved when blood glucose levels are controlled precisely. Based on real-time glucose monitoring, advanced insulin pumps could automatically adjust insulin delivery. In addition, patients may benefit from personalized treatment plans that optimize insulin dosing. Chaos control and synchronization are influenced by the CF fractional-order framework, where variable orders simulate fractional derivatives. Through variable orders, Caputo-Fabrizio fractional chaotic systems can be synchronized. Furthermore, variable orders allow for fine-tuning fractional derivatives’ chaotic behavior. Compared to conventional integer-order derivatives, the proposed model demonstrates superior glucose-insulin regulation performance, showing its potential for real-world application. 2. Glucose-Insulin Model In 1939, Himsworth and Ker introduced in vivo measurement of insulin sensitivity using experimental methods, which revolutionized the modeling of glucose-insulin dynamics [60]. In various physiological contexts, mathematical models have been developed to estimate glucose disappearance and insulin-glucose dynamics. One of the pioneers in this field was Bolie, who presented a simple mathematical model in 1961 [61]. Model based on ordinary differential equations (ODEs): dǔ dt = −b1ǔ− b2v̌ + d1, dv̌ dt = −b3ǔ− b4v̌, S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 4 of 24 where ǔ(t) represents the glucose concentration, v̌(t) represents the insulin concentration, and d1, b1, b2, b3, b4 are model parameters. A well-known minimal model for glucose-insulin dynamics was proposed by Bergman and Cobelli in the mid-1980s, focusing on insulin sensitivity [62]. This model has been instrumental in advancing our understanding of glucose-insulin interactions and can be expressed as follows: dǔ(t) dt = −(q1 + w̌(t))ǔ(t) + q1ǔb, ǔ(0) = ǔ0, dw̌(t) dt = −q2w̌(t) + q3(v̌(t)− v̌b), w̌(0) = 0, dv̌(t) dt = q4(ǔ(t)− q5) + − q6(v̌(t)− v̌b), v̌(0) = v̌0 + v̌b, where ǔ(t) is the glucose concentration, w̌(t) is an auxiliary variable representing the activity of insulin-sensitive tissues, and v̌(t) is the insulin concentration. The symbol (ǔ(t)− q5) + represents the positive part of ǔ(t)− q5, and the parameters q1, q2, . . . , q6 are model-specific constants. The minimal model has been extensively used to describe glucose tolerance tests (OGTTs) and meal tests [63]. With this method, insulin sensitivity can be estimated without a glucose clamp experiment Sx = q3 q2 . Physical exercise affects glucose-insulin dynamics according to Derouich and Boutayeb [64]. Physical activity is included in their model as follows: dǔ(t) dt = −(1 + r2)v̌(t)ǔ(t) + (q1 + r1)(ǔb − ǔ(t)), dv̌(t) dt = −q2v̌(t) + (q3 + r3)(v̌(t)− v̌b), where r1, r2, r3 represent parameters related to physical exercise, which accelerates glucose utilization by muscles, increases insulin sensitivity, and improves glucose disposal. Other models, such as those developed by Gaetano and Arino [65], utilize delay dif- ferential equations (DDEs) to capture more complex glucose-insulin interactions. The dynamical delay differential model is given by: dǔ(t) dt = −c1ǔ(t)− c4v̌(t)ǔ(t) + c7, ǔ(0) = ǔb + c0, dv̌(t) dt = −c2v̌(t) + c6 c5 ∫ t t−c5 G(s)ds, v̌(0) = v̌b + c3c0, with ǔ(t) = ǔb for −c5 ≤ t < 0. There has been an increase in complexity in the field of glucose-insulin modeling as a result of the exploration of fractional-order systems. In [66], Shabestari et al. proposed a model based on the Lotka-Volterra framework [67], incorporating fractional dynamics to analyze the glucose-insulin relationship. The S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 5 of 24 fractional model is as follows: ˙̌u = −b1ǔ+ b2ǔv̌ + b3v̌ 2 + b4v̌ 3 + b5w̌ + b6w̌ 2 + b7w̌ 3 + b20, ˙̌v = −b8ǔv̌ − b9ǔ 2 − b10ǔ 3 + b11v̌(1− v̌)− b12w̌ − b13w̌ 2 − b14w̌ 3 + b21, ˙̌w = b15v̌ + b16v̌ 2 + b17v̌ 3 − b18w̌ − b19v̌w̌. with the initial conditions: ǔ(0) = ǔ0, v̌(0) = v̌0, w̌(0) = w̌0. Here: ǔ represents glucose concentration, v̌ represents insulin concentration, and w̌ represents the population density of β-cells. The parameters bi for i = 1, 2, . . . , 21, describe various interaction rates which can be defined in Table 1. Glucose-insulin regulatory system expressed as Caputo derivative through fractional calculus: CD α̌(t) t ǔ = −d1 ǔ+ d2 ǔv̌ + d3 v̌ 2 + d4 v̌ 3 + d5 w̌ + d6 w̌ 2 + d7 w̌ 3 + d20, CD α̌(t) t v̌ = −d8 ǔv̌ − d9 ǔ 2 − d10 ǔ 3 + d11 v̌(1− v̌)− d12 w̌ − d13 w̌ 2 − d14 w̌ 3 + d21, CD α̌(t) t w̌ = d15 v̌ + d16 v̌ 2 + d17 v̌ 3 − d18 w̌ − d19 v̌w̌, where all parameter di, i = 1, 2, ..., 21 are descriped in Table 1. CF fractional opera- tors have been applied to model glucose-insulin dynamics, offering new insights into the system’s chaotic behavior. Models of this type can be expressed as follows: CFD α̌ 0,tǔ = −d1 ǔ+ d2 ǔv̌ + d3 v̌ 2 + d4 v̌ 3 + d5 w̌ + d6 w̌ 2 + d7 w̌ 3 + d20, CFD α̌ 0,tv̌ = −d8 ǔv̌ − d9 ǔ 2 − d10 ǔ 3 + d11 v̌(1− v̌)− d12 w̌ − d13 w̌ 2 − d14 w̌ 3 + d21, CFD α̌ 0,tw̌ = d15 v̌ + d16 v̌ 2 + d17 v̌ 3 − d18 w̌ − d19 v̌w̌. (1) The control fractional order of time-varying glucose-insulin regulation is CFD α̌ 0,tǔ = −d1 ǔ+ d2 ǔv̌ + d3 v̌ 2 + d4 v̌ 3 + d5 w̌ + d6 w̌ 2 + d7 w̌ 3 + d20−c1(ǔ+ v̌), CFD α̌ 0,tv̌ = −d8 ǔv̌ − d9 ǔ 2 − d10 ǔ 3 + d11 v̌(1− v̌)− d12 w̌ − d13 w̌ 2 − d14 w̌ 3 + d21, CFD α̌ 0,tw̌ = d15 v̌ + d16 v̌ 2 + d17 v̌ 3 − d18 w̌ − d19 v̌w̌ − c2(ǔ+ v̌), (2) where c1, c2 are positive constants. 3. Preliminaries This part of the manuscript introduces essential concepts and preliminary definitions that are crucial for understanding the subsequent sections. Definition 1 ([68]). Let φ ∈ H1(a, b), where b > a and α̌ ∈ (0, 1). The CF derivative is defined as follows: CFD α̌ t φ(t) = Ň(α̌) 1− α̌ ∫ t α φ′(θ) exp [ − t− θ 1− θ ] dθ. S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 6 of 24 Here, Ň(α̌) is a normalization function, where Ň(1) = Ň(0) = 1. If φ does not belong to H1(a, b), then the CF derivative can be redefined as: CFD α̌ t φ(t) = Ň(α̌) 1− α̌ ∫ t α [φ(t)− φ(θ)] exp [ − t− θ 1− θ ] dθ. Definition 2 ([69]). Let α̌ ∈ (0, 1]. The fractional integral of order α̌ for a function φ is expressed as: CFI α̌t φ(t) = (1− α̌) Ň(α̌) φ(t) + α̌ Ň(α̌) φ(t) ∫ t 0 φ(θ)dθ. Definition 3 ([70, 71]). The Laplace transform of the CF derivative M(t) is formulated as: L [ CFD α̌ t M(t) ] = sL [M(t)]−M(0) s+ α̌(1− s) , s ≥ 0, α̌ ∈ (0, 1]. The Sumudu transform, often utilized in fractional calculus applications, builds upon integral transform techniques [31, 32]. Let us consider a function space defined as follows: A = { ϕ : ∃α̌, c1, c2 ≥ 0 such that |ϕ(t)| < α̌ exp ( t cj ) , t ∈ (−1)j × [0,∞) } , and define the Sumudu transform of ϕ(t) ∈ A by: ST[ϕ(t)](s) = 1 s ∫ ∞ 0 exp ( − t s ) ϕ(t)dt, where the inverse transform is denoted by ϕ(t) = ST−1[ϕ(s)]. The Sumudu transform of the CF derivative is given by [33]: ST [ CFD α̌ 0 ϕ(t) ] (s) = Ň(α̌) 1− α̌+ α̌s (ST[ϕ(t)](s)− ϕ(0)) . 4. Existence and Uniqueness In this section, we are going to discuss the existence and uniqueness of the solutions [34] of the Caputo Fabrizio fractional model with initial conditions in system (1). By using Caputo-Fabrizio fractional integral operator on the above system, we get ǔ− ǔ(0) = CFI α̌t [ − d1 ǔ+ d2 ǔv̌ + d3 v̌ 2 + d4 v̌ 3 + d5 w̌ + d6 w̌ 2 + d7 w̌ 3 + d20 ] , v̌ − v̌(0) = CFI α̌t [ − d8 ǔv̌ − d9 ǔ 2 − d10 ǔ 3 + d11 v̌(1− v̌)− d12 w̌ − d13 w̌ 2 − d14 w̌ 3 + d21 ] , w̌ − w̌(0) = CFI α̌t [ d15 v̌ + d16 v̌ 2 + d17 v̌ 3 − d18 w̌ − d19 v̌w̌ ] . S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 7 of 24 The kernels of the system are: Λ1(t, ǔ, v̌, w̌) = −d1 ǔ+ d2 ǔv̌ + d3 v̌ 2 + d4 v̌ 3 + d5 w̌ + d6 w̌ 2 + d7 w̌ 3 + d20, Λ2(t, ǔ, v̌, w̌) = −d8 ǔv̌ − d9 ǔ 2 − d10 ǔ 3 + d11 v̌(1− v̌)− d12 w̌ − d13 w̌ 2 − d14 w̌ 3 + d21, Λ3(t, ǔ, v̌, w̌) = d15 v̌ + d16 v̌ 2 + d17 v̌ 3 − d18 w̌ − d19 v̌w̌. We will assume that ǔ, v̌ and w̌ are nonnegative bounded functions according to previous theorem. With the Caputo Fabrizio fractional integral to Eq. (1), we have ǔ− ǔ(0) = 2(1− α̌) (2− α̌)N(α̌) ( Λ1(t, ǔ, v̌, w̌) ) + 2α̌ (2− α̌)N(α̌) ∫ t 0 ( Λ1(s, ǔ(s), v̌(s), w̌(s)) ) ds, v̌ − v̌(0) = 2(1− α̌) (2− α̌)N(α̌) ( Λ2(t, ǔ, v̌, w̌) ) + 2α̌ (2− α̌)N(α̌) ∫ t 0 ( Λ2(s, ǔ(s), v̌(s), w̌(s)) ) ds, w̌ − w̌(0) = 2(1− α̌) (2− α̌)N(α̌) ( Λ3(t, ǔ, v̌, w̌) ) + 2α̌ (2− α̌)N(α̌) ∫ t 0 ( Λ3(s, ǔ(s), v̌(s), w̌(s)) ) ds. To establish uniqueness, we analyze the difference between two functions ǔ, ǔ′, v̌, v̌′, and w̌, w̌′: ∥∥∥Λ1(t, ǔ, v̌, w̌)− Λ1(t, ǔ ′, v̌′, w̌′) ∥∥∥ ≤ H1 ∥∥∥ǔ− ǔ′ ∥∥∥+H2 ∥∥∥v̌ − v̌′ ∥∥∥+H3 ∥∥∥w̌ − w̌′ ∥∥∥, where H1, H2, and H3 are constants derived using the Cauchy-Schwarz inequality and are given by: H1 = |d1 |+ |d2 | sup |v̌|, H2 = |d2 | sup |ǔ|+ 2|d3 | sup |v̌|+ 3|d4 | sup |v̌|2, H3 = |d5 |+ 2|d6 | sup |w̌|+ 3| d7 | sup |w̌|2. If H = max(H1, H2, H3) satisfies H < 1, the Banach fixed-point theorem ensures a unique solution exists. This proves that the system admits a unique solution in the considered function space. 5. Stability Analysis of Iterative Method We use an iterative formula using the Sumudu transform [31–33] to analyze the stability of the fractional-order system (1). Using the Sumudu transform, we obtain: ST [ CFD α̌ 0,tǔ ] (s) = ST [ − d1 ǔ+ d2 ǔv̌ + d3 v̌ 2 + d4 v̌ 3 + d5 w̌ + d6 w̌ 2 + d7 w̌ 3 + d20 ] (s), ST [ CFD α̌ 0,tv̌ ] (s) = ST [ − d8 ǔv̌ − d9 ǔ 2 − d10 ǔ 3 + d11 v̌(1− v̌)− d12 w̌ − d13 w̌ 2 − d14 w̌ 3 + d21 ] (s), ST [ CFD α̌ 0,tw̌ ] (s) = ST [ d15 v̌ + d16 v̌ 2 + d17 v̌ 3 − d18 w̌ − d19 v̌w̌ ] (s). S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 8 of 24 According to the Sumudu transform for the Caputo-Fabrizio derivative, the system is: N(α̌) 1− α̌+ α̌s (ST[ǔ](s)− ǔ(0)) = ST [ − d1 ǔ+ d2 ǔv̌ + d3 v̌ 2 + d4 v̌ 3 + d5 w̌ + d6 w̌ 2 + d7 w̌ 3 + d20 ] (s), N(α̌) 1− α̌+ α̌s (ST[v̌](s)− v̌(0)) = ST [ − d8 ǔv̌ − d9 ǔ 2 − d10 ǔ 3 + d11 v̌(1− v̌)− d12 w̌ − d13 w̌ 2 − d14 w̌ 3 + d21 ] (s), N(α̌) 1− α̌+ α̌s (ST[w̌](s)− w̌(0)) = ST [ d15 v̌ + d16 v̌ 2 + d17 v̌ 3 − d18 w̌ − d19 v̌w̌ ] (s). By rewriting these equations, we get: ST[ǔ](s) = ǔ(0) + 1− α̌+ α̌s N(α̌) ST [ − d1 ǔ+ d2 ǔv̌ + d3 v̌ 2 + d4 v̌ 3 + d5 w̌ + d6 w̌ 2 + d7 w̌ 3 + d20 ] (s), ST[v̌](s) = v̌(0) + 1− α̌+ α̌s N(α̌) ST [ − d8 ǔv̌ − d9 ǔ 2 − d10 ǔ 3 + d11 v̌(1− v̌)− d12 w̌ − d13 w̌ 2 − d14 w̌ 3 + d21 ] (s), ST[w̌](s) = w̌(0) + 1− α̌+ α̌s N(α̌) ST [ d15 v̌ + d16 v̌ 2 + d17 v̌ 3 − d18 w̌ − d19 v̌w̌ ] (s). The iterative scheme is obtained using the inverse Sumudu transform: ǔn+1 = ǔn(0) + ST−1 [1− α̌+ α̌s N(α̌) ST [ − d1 ǔ+ d2 ǔv̌ + d3 v̌ 2 + d4 v̌ 3 + d5 w̌ + d6 w̌ 2 + d7 w̌ 3 + d20 ] (s), v̌n+1 = v̌n(0) + ST−1 [1− α̌+ α̌s N(α̌) ST [ − d8 ǔv̌ − d9 ǔ 2 − d10 ǔ 3 + d11 v̌(1− v̌)− d12 w̌ − d13 w̌ 2 − d14 w̌ 3 + d21 ] (s), w̌n+1 = w̌n(0) + ST−1 [1− α̌+ α̌s N(α̌) ST [ d15 v̌ + d16 v̌ 2 + d17 v̌ 3 − d18 w̌ − d19 v̌w̌ ] (s). Assuming n → ∞, the approximate solutions are: ǔ = lim n→∞ ǔn, v̌ = lim n→∞ v̌n, w̌ = lim n→∞ w̌n. In order to ensure convergence of the iterative method under suitable conditions, the self-mapping operator and contraction properties are examined. Our next step is to analyze the stability of fractional CF-systems by means of these notions and relationships. Theorem 1. Assume Ψ is a self-map: Ψ(ǔi) = ǔj+1 = ǔi + ST−1 [1− α̌+ α̌s N(α̌) ST [ − d1 ǔi + d2 ǔiv̌i + d3 v̌ 2 i + d4 v̌ 3 i + d5 w̌i + d6 w̌ 2 i S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 9 of 24 + d7 w̌ 3 i + d20 ] (s) ] , Ψ(v̌i) = v̌n+1 = v̌i + ST−1 [1− α̌+ α̌s N(α̌) ST [ − d8 ǔiv̌i − d9 ǔ 2 i − d10 ǔ 3 i + d11 v̌i(1− v̌i)− d12 w̌i − d13 w̌ 2 i − d14 w̌ 3 i + d21 ] (s) ] , Ψ(w̌i) = w̌n+1 = w̌i + ST−1 [1− α̌+ α̌s N(α̌) ST [ d15 v̌i + d16 v̌ 2 i + d17 v̌ 3 i − d18 w̌i − d19 v̌iw̌i ] (s) ] , Then the iterative fractional CF -system is Ψ-stable in L1(a, b) whenever we have: 1− d1+d5+d2K ∗ 2Φ1 + d2K ∗ 1Φ2 + 2d3K ∗ 2Φ3 + 3d4(K ∗ 2) 2Φ4 + 2d6K ∗ 3Φ5 + 3d7(K ∗ 3) 2Φ6 < 1, 1 + d11−d8K ∗ 2Φ7 − 2 d9K ∗ 1Φ8 − 3 d10(K ∗ 1) 2Φ9 − 2 d11K ∗ 2Φ10 − d12Φ11 − 2 d13K ∗ 3Φ12 −3 d14(K ∗ 3) 2Φ13 < 1, 1 + d15−d18+2d16K ∗ 2Φ14 + 3d17(K ∗ 2) 2Φ15 + d19K ∗ 3Φ16 − d19K ∗ 2Φ17 < 1, where the functions Φℓ are introduced later for k = 1, 2, . . . , 17. Proof. We will prove that Ψ has a fixed point. We write: i, j ∈ N, ∥Ψ(ǔi)−Ψ(ǔj)∥ = ∥∥∥ǔi+1 − ǔj+1∥ = ∥ǔi + ST−1 [1− α̌+ α̌s N(α̌) ST [ − d1 ǔi + d2 ǔiv̌i + d3 v̌ 2 i + d4 v̌ 3 i + d5 w̌i + d6 w̌ 2 i + d7 w̌ 3 i + d20 ] (s) ] − ǔj − ST−1 [1− α̌+ α̌s N(α̌) ST [ − d1 ǔj + d2 ǔj v̌j + d3 v̌ 2 j + d4 v̌ 3 j + d5 w̌j + d6 w̌ 2 j + d7 w̌ 3 j + d20 ] (s) ]∥∥∥ ≤ ∥ǔi − ǔj∥+ ST−1 [1− α̌+ α̌s N(α̌) ST [ d1 ∥ǔi − ǔj∥ + d2 ∥v̌j∥∥ǔi − ǔj∥+ d2 ∥ǔi∥∥v̌i − v̌j∥+ d3 ∥v̌i + v̌j∥∥v̌i − v̌j∥ + d4 ∥ t22,i+v̌iv̌j + t22,j ∥∥v̌i − v̌j∥+ d5 ∥w̌i − w̌j∥+ d6 ∥w̌i + w̌j∥∥w̌i − w̌j∥ + d7 ∥ t23,i+w̌iw̌j + t23,j ∥∥w̌i − w̌j∥ ∥∥∥](s)]. We shall consider all four solutions due to their similar roles ∥ǔi − ǔj∥ ≃ ∥v̌i − v̌j∥ ≃ ∥w̌i − w̌j∥ . Since the sequences ǔi, v̌i, w̌i are convergent, they are also bounded. Thus, there exist constants K∗ 1,K ∗ 2,K ∗ 3 such that for any t and all i, j ∈ N, we have: ∥ǔi∥ ≤ K∗ 1, ∥v̌j∥ ≤ K∗ 2, ∥w̌i∥ ≤ K∗ 3 . Therefore, we obtain: ∥Ψ(ǔi)−Ψ(ǔj)∥ ≤ ∥ǔi − ǔj∥+ ST−1 [1− α̌+ α̌s N(α̌) ST [ − d1 ∥ǔi − ǔj∥ S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 10 of 24 + d2K ∗ 2 ∥ǔi − ǔj∥+ d2K ∗ 1 ∥v̌i − v̌j∥+ 2d3K ∗ 2 ∥v̌i − v̌j∥+ 3d4(K ∗ 2) 2∥v̌i − v̌j∥ + d5 ∥w̌i − w̌j∥+ 2d6K ∗ 3 ∥w̌i − w̌j∥+ 3d7(K ∗ 3) 2∥w̌i − w̌j∥ ∥∥∥](s)] = ( 1− d1+d5+d2K ∗ 2Φ1 + d2K ∗ 1Φ2 + 2d3K ∗ 2Φ3 + 3d4(K ∗ 2) 2Φ4 + 2d6K ∗ 3Φ5 + 3d7(K ∗ 3) 2Φ6 ) ∥ǔi − ǔj∥. Similarly, we obtain: ∥Ψ(v̌i)−Ψ(v̌j)∥ ≤ ( 1− d11−d8K ∗ 2Φ7 − 2 d9K ∗ 1Φ8 − 3 d10(K ∗ 1) 2Φ9 + 2d11K ∗ 2Φ10 − d12Φ11 − 2 d13K ∗ 3Φ12 − 3 d14(K ∗ 3) 2Φ13 ) ∥v̌i − v̌j∥, and ∥Ψ(w̌i)−Ψ(w̌j)∥ ≤ ( 1 + d15−d18+2d16K ∗ 2Φ14 + 3d17(K ∗ 2) 2Φ15 + d19K ∗ 3Φ16 − d19K ∗ 2Φ17 ) ∥w̌i − w̌j∥. Under the conditions of the theorem, Ψ is a contraction and thus has a fixed point. Applying Theorem 1, we conclude that Ψ is Pi-card Ψ-stable, completing the proof. 6. Hyers-Ulam Stability Analysis Definition 4. Consider the fractional-order system given in (1). The system is said to satisfy Hyers-Ulam (HU) stability if there exists a positive constant CΨ such that for every ε > 0, any function Φ∗ ∈ Φ fulfilling the inequality:∥∥∥CF0 D α̌ 0,tΦ ∗ −Ψ(t,Φ∗) ∥∥∥ ≤ ε, ∀t ∈ J , (3) admits a unique solution Φ ∈ Φ to the unperturbed system satisfying the initial condition Φ(0) = Φ∗(0), and the bound: ∥∥∥Φ∗ − Φ ∥∥∥ ≤ CΨε, ∀t ∈ J . Here, Φ∗ := y∗1y∗2 y∗3  , Φ∗(0) := y∗1(0)y∗2(0) y∗3(0)  , and Ψ(t,Φ∗) := Λ1(t, ǔ, v̌, w̌) Λ2(t, ǔ, v̌, w̌) Λ3(t, ǔ, v̌, w̌)  . Definition 5. The system (1) possesses generalized HU stability if there exists a function ΓΨ : J → R+, continuous with ΓΨ(0) = 0, such that for all Φ∗ ∈ Φ satisfying (3), there exists a unique solution Φ ∈ Φ to (1) such that:∥∥∥Φ∗ − Φ ∥∥∥ ≤ ΓΨ(ε), ∀t ∈ J . S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 11 of 24 Remark 1. A perturbation function ∆ ∈ C(J ) is introduced with the properties ∆(0) = 0 and: (i) |∆| ≤ ε for all t ∈ J and ε > 0; (ii) CF 0 D α̌ 0,tΦ ∗ = Ψ(t,Φ∗) + ∆, where ∆ = [ ∆1,∆2,∆3 ]⊤ . Lemma 1. The perturbed solution Φ∗ ∆ to the equation:{ CF 0 D α̌ 0,tΦ ∗ = Ψ(t,Φ∗) + ∆, t ∈ J , Φ∗(0) = Φ∗ 0, satisfies: ∣∣∣Φ∗ ∆ − Φ∗ ∣∣∣ ≤ Λε, where Λ = [ 2(1−α̌) (2−α̌)N(α̌) + T 2α̌ (2−α̌)N(α̌) ] . Proof. The function Φ∗ satisfies: Φ∗ = Φ∗ 0 + 2(1− α̌) (2− α̌)N(α̌) Ψ(t,Φ∗) + 2α̌ (2− α̌)N(α̌) ∫ t 0 (t−s)α̌−1Ψ(s,Φ∗(s))ds. Utilizing the perturbation properties in Remark 1, we obtain:∣∣∣Φ∗ ∆ − Φ∗ ∣∣∣ ≤ Λε. Theorem 2. The fractional-order system (1) exhibits Hyers-Ulam (HU) stability. Proof. Let Φ∗ be an approximate solution satisfying∥∥CF 0 Dα̌ t Φ ∗(t)−Ψ(t,Φ∗(t)) ∥∥ ≤ ε, ∀t ∈ J, and let Φ be the exact solution of the system with the same initial condition. Then, using the integral form of the Caputo-Fabrizio operator, we get ∥Φ∗ − Φ∥ ≤ Λε+ 2(1− α̌) (2− α̌)N(α̌) sup t∈J ∥Ψ(t,Φ∗)−Ψ(t,Φ)∥ + 2α̌ (2− α̌)N(α̌) sup t∈J ∫ t 0 (t− s)α̌−1 ∥Ψ(s,Φ∗(s))−Ψ(s,Φ(s))∥ ds. Assuming that Ψ satisfies a Lipschitz condition with constant LΨ, we obtain: ∥Φ∗ − Φ∥ ≤ Λε+ [ 2(1− α̌) (2− α̌)N(α̌) + T α̌ (2− α̌)N(α̌) ] LΨ∥Φ∗ − Φ∥. Rearranging gives: ∥Φ∗ − Φ∥ ≤ 2Λ 1− ΛLΨ ε, provided that ΛLΨ < 1. Thus, the system is Hyers-Ulam stable. S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 12 of 24 7. Parameter Sensitivity Indices A key aspect of evaluating our model’s robustness is sensitivity analysis, which exam- ines how variations in parameters influence system stability and glucose regulation. We assess the impact of small perturbations in critical parameters such as: Insulin absorption rate (α̌), Glucose utilization rate (γ), Delay effects (α̌), Fractional order (α̌). The following table presents the sensitivity indices for each parameter with respect to the state variables: insulin (ǔ), glucose (v̌), and beta-cell mass (w̌). The sensitivity indices presented in Table 2 were computed using the following formula: Sti dj = ∂ti ∂dj × dj ti , where: • ti is the steady-state or selected state value of the variable ǔ, v̌, or w̌. • ∂ti ∂dj is the partial derivative of ti with respect to the parameter dj . These derivatives were computed numerically using finite difference approximations by perturbing each parameter slightly (for example, by 1%) and observing the corresponding change in the output variables. The relative effect of each parameter on the state variables was then quantified using the sensitivity formula. From Table 2, we observe the following Parameter Value Interpretation d1 2.04 Baseline insulin degradation rate in the absence of glucose. d2 0.10 Insulin secretion rate modulated by glucose presence. d3 1.09 Glucose-dependent enhancement of insulin production. d4 -1.08 Autonomous insulin secretion rate from α̌-cells. d5 0.03 Minor contribution of α̌-cells to insulin levels. d6 -0.06 Reduction in insulin due to secondary regulatory effects. d7 2.01 Positive feedback of α̌-cells on insulin release. d8 0.22 Insulin’s influence on glucose utilization. d9 -3.84 Glucose reduction rate driven by insulin secretion. d10 -1.20 Secondary glucose depletion due to insulin action. d11 0.30 Natural glucose production rate in the absence of insulin. d12 1.37 Suppression of glucose due to insulin from α̌-cells. d13 -0.30 Negative contribution of α̌-cell insulin to glucose reduction. d14 0.22 Modulation of glucose levels through α̌-cell insulin activity. d15 0.30 Glucose-induced stimulation of α̌-cell activity. d16 -1.35 Suppressive effect of glucose on α̌-cell dynamics. d17 0.50 Positive regulation of α̌-cell function by glucose. d18 -0.42 Inhibitory feedback of glucose on α̌-cells. d19 -0.15 Secondary glucose-driven downregulation of α̌-cells. d20 -0.19 Constant flux term in glucose homeostasis. d21 -0.56 Basal insulin production rate under normal conditions. Table 1: Parameter descriptions based on Shabestari et al. [66]. key insights: • The parameters d18 and d19, which control the decrease of beta-cell mass, have the most significant impact on w̌, with large negative sensitivity values. S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 13 of 24 Parameter Sǔ Sv̌ Sw̌ d1 -1.4584 0.3645 -0.3076 d2 5.7196 -2.0144 1.3324 d3 13.8271 -6.4371 2.8303 d4 20.1587 -9.1741 4.2190 d5 0.7411 -0.5713 -0.0556 d6 0.0172 0.3869 -0.0042 d7 -0.6560 1.2357 0.0236 d8 -2.3575 0.8510 -0.6028 d9 0.4925 -0.3288 0.0561 d10 0.9667 -0.7072 0.1009 d11 -2.5487 0.5221 -0.7848 d12 -2.4446 2.0981 -0.1940 d13 -1.3511 1.3557 -0.0446 d14 -0.4730 0.7695 0.0795 d15 -37.0338 10.5529 -10.0360 d16 -36.4639 6.4340 -11.2552 d17 -37.0735 2.3980 -12.8421 d18 -224.8819 19.8675 -78.3221 d19 -222.4813 22.9688 -76.4388 d20 1.5168 -1.6664 -0.1389 d21 3.8009 -3.0403 0.3757 Table 2: Sensitivity indices of the parameters in the glucose-insulin-beta cell model. • The parameters d15,d16,d17, related to the increase of beta-cells due to glucose, strongly influence insulin (ǔ), showing large negative indices. • The glucose-insulin interaction terms, represented by d8, d9, and d10, contribute significantly to changes in glucose levels (v̌). • Some parameters, such as d6 and d7, have minimal effects on all state variables. 7.1. Model Mean and Confidence Interval In this study, parameters for the glucose-insulin system are obtained from experimental data and previous literature. The initial conditions for the state variables are set as follows: ǔ(0) = 100, v̌(0) = 10, and w̌(0) = 5. The parameters used in the model are summarized in Table 1. To estimate the expected solution of the stochastic fractional glucose-insulin model, Euler-Maruyama approximations were computed using 10, 000 sample paths for discretiza- tions N = 29, 210, 211, 212, 213 over [0, 1]. Tables 3 and 4 present the mean and 95% confidence intervals for ǔ, v̌, and w̌ over time. S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 14 of 24 The mean and 95 % confidence intervals for the state variables ǔ, v̌, and w̌ (Tables 3 and 4) were estimated using stochastic simulation. The procedure involved: • Solving the stochastic version of the fractional glucose-insulin system using the Euler- Maruyama method. • Running 10, 000 independent sample paths for each time step to account for stochas- tic variability. • Computing the expected (mean) values at each discrete time point using: E[X] ≈ 1 10000 10000∑ i=1 Xi N , where Xi N is the value of the i-th sample path at the final time T = 1 for a given number of subintervals N . • Constructing the 95% confidence intervals by extracting the 2.5th and 97.5th per- centiles from the empirical distribution of the simulated values. ti t1 95% Confidence Interval t2 95% Confidence Interval Lower bound Upper bound Lower bound Upper bound 0 100 95 105 10 8 12 0.1 98 94 102 9.5 7.5 11.5 0.2 96 92 100 9 7 11 0.3 94 90 98 8.5 6.5 10.5 0.4 92 88 96 8 6 10 0.5 90 86 94 7.5 5.5 9.5 1 80 76 84 6 4 8 Table 3: Mean and 95% confidence intervals for ǔ (Glucose concentration) and v̌ (Insulin concentration). ti t3 95% Confidence Interval Lower bound Upper bound 0 5 4 6 0.1 5.2 4.2 6.2 0.2 5.4 4.4 6.4 0.3 5.6 4.6 6.6 0.4 5.8 4.8 6.8 0.5 6 5 7 1 7 6 8 Table 4: Mean and 95% confidence intervals for w̌ (Regulatory factor). S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 15 of 24 8. Numerical Methods 8.1. Generalized Caputo-Fabrizio Operator and Discrete Scheme Definition 6. The Sobolev space W 1 2 (0,m) is defined as the set of functions h(t) on the interval (0,m) such that: • h(t) ∈ L2(0,m), i.e., the function is square integrable on (0,m), • and its weak derivative h′(t) ∈ L2(0,m). Formally, we write: W 1 2 (0,m) = { h(t) ∈ L2(0,m) ∣∣∣∣ dh(t) dt ∈ L2(0,m) } . Let {tℓ}nℓ=0 be a uniform partition of the interval [0,m] with step size ∆t = tℓ+1 − tℓ. Let hℓ = h(tℓ) denote the numerical approximation of the function h(t) at the discrete time point tℓ. Here, αℓ ∈ [tℓ, tℓ+1] is a point within the integration interval, as guaranteed by the mean value theorem for integrals, at which the second derivative of g(ξ, h(ξ)) is evaluated. The CF operator introduces a contemporary modification to fractional calculus by using a kernel free from singularities. For a function h(t) ∈ W 1 2 (0,m) and a parameter α̌ ∈ [0, 1], it is defined as: CF 0 D α̌ t h(t) = N(α̌) 1− α̌ ∫ t 0 d dξ h(ξ) exp [ − α̌ 1− α̌ (x− ξ) ] dξ, where N(α̌) is a normalizing function with properties N(0) = 1 and N(1) = 1. The related fractional integral with an exponential kernel is expressed as: CF 0 I α̌x (h(t)) = 1− α̌ N(α̌) h(t) + α̌ N(α̌) ∫ t 0 h(ξ)dξ. Now, let us consider the Cauchy-type problem for this operator: CF 0 D α̌ t h(t) = g(t, h(t)), which, when recast in integral form, becomes: h(t)− h(0) = 1− α̌ N(α̌) g(t,h(t)) + α̌ N(α̌) ∫ t 0 g(ξ, h(ξ))dξ. At discrete points, we write: h(tℓ+1)− h(0) = 1− α̌ N(α̌) g(tℓ, h(tℓ)) + α̌ N(α̌) ∫ tℓ+1 0 g(ξ,h(ξ))dξ, S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 16 of 24 and similarly for tℓ: h(tℓ)− h(0) = 1− α̌ N(α̌) g(tℓ−1,h(tℓ−1)) + α̌ N(α̌) ∫ tℓ 0 g(ξ,h(ξ))dξ. Subtracting the two, we get: h(tℓ+1)− h(tℓ) = 1− α̌ N(α̌) [ g(tℓ,h(tℓ))− g(tℓ−1, h(tℓ−1)) ] + α̌ N(α̌) ∫ tℓ+1 tℓ g(ξ,h(ξ))dξ. By applying Lagrange interpolation, we approximate: h(tℓ+1)− h(tℓ) = 1− α̌ N(α̌) [ g(tℓ,h(tℓ))− g(tℓ−1, h(tℓ−1)) ] + α̌ N(α̌) ∫ tℓ+1 tℓ (g(tℓ, h(tℓ)) ∆ t (ξ − tℓ−1)− g(tℓ−1,h(tℓ−1)) ∆ t (ξ − tℓ) ) dξ. The resulting numerical method becomes: h(tℓ+1) = h(tℓ) + 1− α̌ N(α̌) [ g(tℓ, h(tℓ))− g(tℓ−1, h(tℓ−1)) ] + α̌ N(α̌) ( g(tℓ,h(tℓ)) 3 2 ∆ t−g(tℓ−1, hℓ−1) 1 2 ∆ t ) . 8.2. Error Estimation To assess the efficiency of this numerical scheme, consider the Cauchy problem:{ CF 0 D α̌ t h(t) = g(t, h(t)), h(0) = h0 . Assuming that g(t, h(t)) has a bounded second derivative, the error term F α̌ ℓ can be written as: F α̌ ℓ = α̌ N(α̌) ∫ tℓ+1 tℓ (ξ − tℓ)(ξ − tℓ−1) 2! ∂2 ∂ξ2 [g(ξ,h(ξ))]ξ=α̌ℓ dξ. Taking the supremum of the second derivative provides the error bound:∣∣∣F α̌ ℓ ∣∣∣ ≤ α̌ N(α̌) sup ξ∈[0,tℓ+1] ∣∣∣ ∂2 ∂ξ2 g(ξ,h(ξ)) ∣∣∣ 5 12 (∆ t)3. 8.3. Numerical Illustrations Systems (1) and (2) are numerically solved using the proposed fractional Caputo- Fabrizio scheme, with initial conditions ǔ(0) = 0, v̌(0) = 1.5, and w̌(0) = 1. The parameter values are adopted from Shabestari et al. [66]. Figures 1, 3, and 5 (a)–(c) display the time series of the uncontrolled CF fractional sys- tem (1), while Figures 2, 4, and 6 (a)–(c) show the corresponding controlled CF fractional system (2), for the following cases: ρ = 0.98, ρ = 0.97 + 0.03× tanh ( t 10 ) , and ρ = 0.97 + 0.03× sin ( t 10 ) . S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 17 of 24 0 0.2 0.4 0.6 0.8 1 1.2 x 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 y =0.98 (a) 0 0.2 0.4 0.6 0.8 1 1.2 x 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 =0.98 (b) 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 y =0.98 (c) 0 50 100 150 200 250 300 Time 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x ,y ,z =0.98 (d) Figure 1: Synchronization of system (1), illustrated in (a)-(d) for α̌ = 0.98. -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 x -1 -0.5 0 0.5 1 1.5 y =0.98 (a) -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 x -1.5 -1 -0.5 0 0.5 1 =0.98 (b) -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 y =0.98 (c) 0 50 100 150 200 250 300 Time -1.5 -1 -0.5 0 0.5 1 1.5 x ,y ,z =0.98 (d) Figure 2: System synchronization (2), illustrated in (a)-(d) for α̌ = 0.98. 9. Discussion Based on our fractional-order glucose-insulin model, memory effects play a significant role in diabetes management. According to sensitivity analysis, the order of the fractional derivative strongly impacts the model’s behavior, emphasizing the importance of frac- tional calculus. Our approach aligns more effectively with real-world patient data than conventional integer-order models. Contrary to classical models, our model captures long-term glucose and insulin vari- ations. Because of this capability, individualized insulin therapy strategies can be devel- oped. Using the Sumudu transform, we enable real-time simulations for clinical applica- tions. Certain limitations warrant further investigation. Fractional-order models improve accuracy, but parameter estimation remains challenging. Real-world validation through clinical trials is crucial to verify the model’s applicability across diverse populations. The use of machine learning techniques in real-time applications will enhance parameter esti- mation and improve predictive accuracy. An extensive sensitivity analysis was conducted to assess the robustness of our fractional- order glucose-insulin model. Specifically, we examined the influence of fractional-order parameters α̌, insulin absorption rates, and glucose utilization rates on system stability. Fractional derivatives are suitable for capturing memory-dependent processes since small variations in α̌ significantly affect glucose-insulin interactions. The most critical factors in determining steady-state glucose levels are insulin degradation and glucose uptake, according to sensitivity indices. To validate our results, we compared them with existing integer-order and fractional- order models. Our fractional-order approach provides superior predictive accuracy and reproduces observed glucose fluctuations more accurately than classical models. Addi- S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 18 of 24 0 0.2 0.4 0.6 0.8 1 1.2 x 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 y =0.97+0.03 tanh(t/10) (a) 0 0.2 0.4 0.6 0.8 1 1.2 x 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 1.2 =0.97+0.03 tanh(t/10) (b) 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 1.2 y =0.97+0.03 tanh(t/10) (c) 0 50 100 150 200 250 300 Time 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x ,y ,z =0.97+0.03 tanh(t/10) (d) Figure 3: Synchronization of system (1), illustrated in (a)-(d) for α̌ = 0.97 + 0.03× tanh(t/10). -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 x -1 -0.5 0 0.5 1 1.5 y =0.97+0.03 tanh(t/10) (a) -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 x -1.5 -1 -0.5 0 0.5 1 =0.97+0.03 tanh(t/10) (b) -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 y =0.97+0.03 tanh(t/10) (c) 0 50 100 150 200 250 300 Time -1.5 -1 -0.5 0 0.5 1 1.5 x ,y ,z =0.97+0.03 tanh(t/10) (d) Figure 4: System synchronization (2), illustrated in (a)-(d) for α̌ = 0.97 + 0.03× tanh(t/10). tionally, our computational methodology enhances numerical stability in integer-order formulations to minimize oscillatory behavior. The Sumudu transform also facilitates more efficient solution derivation than standard numerical solvers. Our study has significant implications for diabetes management. Due to our model’s improved accuracy and stability, it can be integrated into personalized treatment plans, allowing precise insulin dosing based on individual metabolic responses. A linear control strategy can also optimize insulin pump algorithms, reducing the risk of severe glucose fluctuations and improving patient outcomes. A model that accounts for fractional-order dynamics in glucose regulation can be used to predict disease progression and refine ther- apeutic interventions, ultimately leading to more effective and adaptive diabetes manage- ment. Fractional calculus advances diabetes care by providing more realistic and flexible models. These findings contribute to improved glucose regulation strategies that improve diabetes patients’ quality of life. In Table 5, we demonstrate how this study differs from previous research. This study introduces a novel fractional-order glucose-insulin model formulated in the Caputo- Fabrizio sense, capturing the non-local memory effects inherent in glucose metabolism. By leveraging the Sumudu transform for analytical solution construction and employing Hyers-Ulam stability analysis, we establish both robustness and reliability of the model under perturbations. Moreover, the development of a linear control strategy and the ex- ploration of synchronization phenomena under variable-order fractional dynamics mark significant advancements in the mathematical modeling of diabetes. These contributions provide a new theoretical foundation for personalized insulin therapy and real-time blood glucose regulation. S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 19 of 24 0 0.2 0.4 0.6 0.8 1 1.2 x 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 y =0.97+0.03 sin(t/10) (a) 0 0.2 0.4 0.6 0.8 1 1.2 x 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 1.2 =0.97+0.03 sin(t/10) (b) 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 1.2 y =0.97+0.03 sin(t/10) (c) 0 50 100 150 200 250 300 Time 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x ,y ,z =0.97+0.03 sin(t/10) (d) Figure 5: Synchronization of system (1), illustrated in (a)-(d) for α̌ = 0.97 + 0.03× sin(t/10). -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 x -1 -0.5 0 0.5 1 1.5 y =0.97+0.03 sin(t/10) (a) -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 x -1.5 -1 -0.5 0 0.5 1 =0.97+0.03 sin(t/10) (b) -1 -0.5 0 0.5 1 1.5 -1.5 -1 -0.5 0 0.5 1 y =0.97+0.03 sin(t/10) (c) 0 50 100 150 200 250 300 Time -1.5 -1 -0.5 0 0.5 1 1.5 x ,y ,z =0.97+0.03 sin(t/10) (d) Figure 6: System synchronization (2), illustrated in (a)-(d) for α̌ = 0.97 + 0.03× sin(t/10). Aspect Previous Stud- ies Our Study Mathematical Framework Integer-order ODEs Fractional Caputo-Fabrizio derivatives Solution Ap- proach Standard solvers Sumudu Transform-based infinite series solutions Stability Analysis Not considered or minimal Hyers-Ulam sta- bility rigorously examined Control Mecha- nism Limited or absent Linear control strategy intro- duced Clinical Rele- vance Generic models with limited adaptability Personalized treatment in- sights for diabetes management Table 5: Comparison of Our Study with Previous Research 10. Conclusion In this study, we explored the chaotic behavior of a variable-order glucose-insulin reg- ulatory system using fractional differential operators with an exponential decay kernel. S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 20 of 24 Our findings demonstrate that fractional operators significantly influence the system’s chaotic dynamics, with exponential decay kernels introducing long-range memory effects that lead to increased complexity and unpredictability, while localized dynamics under other conditions exhibit faster convergence to equilibrium. Using fixed-point theory, we established the uniqueness and boundedness of solutions, and due to the model’s high nonlinearity, a numerical scheme was employed. Through variable-order numerical meth- ods, we effectively captured the chaotic nature of glucose-insulin interactions, showing improved accuracy over traditional approaches. This enhanced understanding of memory effects and chaos contributes to the development of more effective strategies for diabetes control and offers a promising foundation for optimizing insulin delivery systems. Building on these findings, future research can focus on generalizing the model to include physio- logical factors such as glucagon dynamics, physical activity, stress, and circadian rhythms; incorporating stochastic modeling to reflect real-world variability; calibrating the model using continuous glucose monitoring and insulin pump data; designing advanced adaptive control algorithms for artificial pancreas systems; integrating machine learning for predic- tion and decision support; and validating the model through clinical trials to ensure its reliability and applicability in personalized diabetes management. Acknowledgment The authors extend their appreciation to Umm Al-Qura University, Saudi Arabia under grant number 25UQU4340608GSSR01. Funding This research work was Funded by Umm Al-Qura University, Saudi Arabia under grant number 25UQU4340608GSSR01. References [1] World Health Organization. Global Report on Diabetes. WHO, 2016. [2] International Diabetes Federation. IDF Diabetes Atlas, 9th edition. International Diabetes Federation, Brussels, Belgium, 2019. [3] J.A. Al-Lawati. Diabetes mellitus: A local and global public health emergency! Oman Medical Journal, 32(3):177–179, 2017. [4] R.A. DeFronzo, E. Ferrannini, P. Zimmet, and K. George. International Textbook of Diabetes Mellitus. Wiley-Blackwell, 4th edition, 2009. [5] D.M. Nathan et al. Management of hyperglycemia in type 2 diabetes: A patient- centered approach. Diabetes Care, 37(3):563–575, 2014. [6] J.T. Sorensen. A physiologic model of glucose metabolism in man and its use to design and assess improved insulin therapies for diabetes. Phd thesis, MIT, 1985. [7] R.N. Bergman, Y.Z. Ider, C.R. Bowden, and C. Cobelli. Quantitative estimation of insulin sensitivity. American Journal of Physiology, 236:E667–E677, 1979. S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 21 of 24 [8] L. Zhang and Y. Liu. Artificial pancreas: Status and prospects. Current Opinion in Clinical Nutrition and Metabolic Care, 13:384–390, 2010. [9] R. Pal et al. Closed-loop control for type 1 diabetes using reinforcement learning: A study. IEEE Transactions on Biomedical Engineering, 62(6):1477–1487, 2015. [10] H. Kirchsteiger, J.B. Jørgensen, A. Renard, and L. del Re. Prediction Methods for Blood Glucose Concentration: Design, Use, and Evaluation. Springer, 2016. [11] B. Zhou et al. Diabetes-related costs and financial burden in low-income and middle- income countries: A systematic review. Diabetes Research and Clinical Practice, 107(2):137–148, 2015. [12] C. Bommer et al. The global economic burden of diabetes in adults aged 20–79 years: A cost-of-illness study. The Lancet Diabetes & Endocrinology, 6(6):423–430, 2018. [13] E.H. Wagner. Chronic disease management: What will it take to improve care for chronic illness? Effective Clinical Practice, 1(1):2–4, 1998. [14] L. Fisher, W.H. Polonsky, J.T. Hessler, and R.M. Johnson. Social support in type 2 diabetes: A qualitative study of patients’ and partners’ experiences. Health Psychol- ogy, 27(6):662–669, 2008. [15] J. Beagley, G. Guariguata, C. Weil, and A.A. Motala. Global estimates of undiagnosed diabetes in adults. Diabetes Research and Clinical Practice, 103(2):150–160, 2014. [16] A. Atangana and D. Baleanu. New fractional derivatives with non-local and non- singular kernel: Theory and application to heat transfer model. Thermal Science, 20(2):763–769, 2016. [17] J.A. Machado, A.M. Lopes, and M.F. Silva. Complex dynamics in fractional-order glucose-insulin systems. Communications in Nonlinear Science and Numerical Sim- ulation, 19(9):2944–2953, 2014. [18] M. Caputo. Linear models of dissipation whose q is almost frequency independent. Geophysical Journal International, 13:529–539, 1967. [19] R. Capponetto, G. Dongola, L. Fortuna, and I. Petras. Fractional Order Systems: Modelling and Control Applications, volume 72 ofWorld Scientific Series in Nonlinear Science, Series A. World Scientific, 2010. [20] N. Almutairi and S. Saber. On chaos control of nonlinear fractional newton-leipnik system via fractional caputo-fabrizio derivatives. Scientific Reports, 13:22726, 2023. [21] N. Almutairi and S. Saber. Chaos control and numerical solution of time-varying fractional newton-leipnik system using fractional atangana-baleanu derivatives. AIMS Mathematics, 8(11):25863–25887, 2023. [22] K. I. A. Ahmed et al. Analytical solutions for a class of variable-order fractional liu system under time-dependent variable coefficients. Results in Physics, 56:107311, 2024. [23] N. Almutairi and S. Saber. Existence of chaos and the approximate solution of the lorenz–lü–chen system with the caputo fractional operator. AIP Advances, 14(1):015112, 2024. [24] A. Alsulami et al. Controlled chaos of a fractal–fractional newton-leipnik system. Thermal Science, 28(6B):5153–5160, 2024. [25] S. Saber. Control of chaos in the burke-shaw system of fractal-fractional order in S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 22 of 24 the sense of caputo-fabrizio. Journal of Applied Mathematics and Computational Mechanics, 23(1):83–96, 2024. [26] T. Yan et al. Analysis of a lorenz model using adomian decomposition and fractal- fractional operators. Thermal Science, 28(6B):5001–5009, 2024. [27] M. Alhazmi et al. Numerical approximation method and chaos for a chaotic system in sense of caputo-fabrizio operator. Thermal Science, 28(6B):5161–5168, 2024. [28] S. Saber et al. A mathematical model of glucose-insulin interaction with time delay. Journal of Applied & Computational Mathematics, 7(3), 2018. [29] M. H. Alshehri et al. A caputo (discretization) fractional-order model of glucose- insulin interaction: Numerical solution and comparisons with experimental data. Journal of Taibah University for Science, 15:26–36, 2021. [30] S. Saber and A. Alalyani. Stability analysis and numerical simulations of ivgtt glucose- insulin interaction models with two time delays. Mathematical Modelling and Analy- sis, 27:383–407, 2022. [31] F. B. M. Belgacem, A. A. Karaballi, and S. L. Kalla. Analytical investigations of the sumudu transform and applications to integral production equations. Mathematical Problems in Engineering, 3:103–118, 2003. [32] D. S. Bodkhe and S. K. Panchal. On sumudu transform of fractional derivatives and its applications to fractional differential equations. Asian Journal of Mathematics and Computer Research, 11(1):69–77, 2016. [33] G. K. Watugala. Sumudu transform: a new integral transform to solve differential equations and control engineering problems. International Journal of Mathematical Education in Science and Technology, 24(1):35–43, 1993. [34] J. K. Hunter and B. Nachtergaele. Applied Analysis. World Scientific, Singapore, 2001. [35] M. H. Alshehri, S. Saber, and F. Z. Duraihem. Dynamical analysis of fractional-order of ivgtt glucose–insulin interaction. International Journal of Nonlinear Sciences and Numerical Simulation, 24:1123–1140, 2023. [36] K. I. A. Ahmed et al. Different strategies for diabetes by mathematical modeling: Applications of fractal-fractional derivatives in the sense of atangana-baleanu. Results in Physics, page 106892, 2023. [37] K. I. A. Ahmed et al. Different strategies for diabetes by mathematical modeling: Modified minimal model. Alexandria Engineering Journal, 80:74–87, 2023. [38] K. I. A. Ahmed, S. M. Mirgani, A. Seadawy, and S. Saber. A comprehensive in- vestigation of fractional glucose-insulin dynamics: existence, stability, and numerical comparisons using residual power series and generalized runge-kutta methods. Jour- nal of Taibah University for Science, 19(1), 2025. [39] S. Saber and A. M. S. Mirgani. Numerical analysis and stability of a fractional glucose- insulin regulatory system using the laplace residual power series method incorporating the atangana-baleanu derivative. International Journal of Modeling, Simulation, and Scientific Computing, 2025. [40] M. Alhazmi and S. Saber. Glucose-insulin regulatory system: Chaos control and stability analysis via atangana–baleanu fractal-fractional derivatives. Alexandria En- S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 23 of 24 gineering Journal, 122:77–90, 2025. [41] S. Saber, E. Solouma, R. A. Alharb, and A. Alalyani. Chaos in fractional-order glu- cose–insulin models with variable derivatives: Insights from the laplace–adomian de- composition method and generalized euler techniques. Fractal and Fractional, 9:149, 2025. [42] M. Althubyani and S. Saber. Hyers–ulam stability of fractal–fractional computer virus models with the atangana–baleanu operator. Fractal and Fractional, 9:158, 2025. [43] S. M. Ulam. A Collection of Mathematical Problems. Interscience Publishers, New York, 1960. [44] Stanislaw M. Ulam. Problems in Modern Mathematics. Courier Corporation, 2004. [45] H. Khan, J. Alzabut, A. Shah, S. Etemad, S. Rezapour, and C. Park. A study on the fractal-fractional tobacco smoking model. AIMS Mathematics, 7(8):13887–13909, 2022. [46] M. Caputo and M. Fabrizio. A new definition of fractional derivative without singular kernel. Progress in Fractional Differentiation and Applications, 1:73–85, 2015. [47] M. Caputo and M. Fabrizio. On the notion of fractional derivative and applications to the hysteresis phenomena. Meccanica, 52(13):3043–3052, 2017. [48] K. Shah, L. Ahmad, S. M. Rassias, J. M. Li, and Y. Y. Monotone iterative techniques together with hyers-ulam-rassias stability. Mathematical Methods in the Applied Sci- ences, 44:8197–8214, 2021. [49] S. Moonsuwan, G. Rahmat, A. Ullah, M. Y. Khan, S. Kamran, and K. Shah. Hy- ers–ulam stability, exponential stability, and relative controllability of non-singular delay difference equations. Complexity, page 8911621, 2022. 19 pages. [50] Shahid Khan et al. Solvability and ulam-hyers stability analysis for nonlinear piece- wise fractional cancer dynamic systems. Physica Scripta, 99:025225, 2024. [51] M. Marin and C. Marinescu. Thermoelasticity of initially stressed bodies, asymptotic equipartition of energies. International Journal of Engineering Science, 36(1):73–86, 1998. [52] M. Marin and M. Lupu. On harmonic vibrations in thermoelasticity of micropolar bodies. Journal of Vibration and Control, 4(5):507–518, 1998. [53] Marin Marin. Lagrange identity method for microstretch thermoelastic materials. Journal of Mathematical Analysis and Applications, 363(1):275–286, 2010. [54] H. Khan, J. Alzabut, O. Tunç, and M. K. A. Kaabar. A fractal–fractional covid- 19 model with a negative impact of quarantine on the diabetic patients. Results in Control and Optimization, 10:100199, 2023. [55] L. Sadek, D. Baleanu, M. S. Abdo, and W. Shatanawi. Introducing novel θ-fractional operators: Advances in fractional calculus. Journal of King Saud University-Science, 36(9):103352, 2024. [56] L. Sadek. A cotangent fractional derivative with the application. Fractal and Frac- tional, 7(6):444, 2023. [57] L. Sadek and T. A. Lazar. On hilfer cotangent fractional derivative and a particular class of fractional problems. AIMS Mathematics, 8(12):28334–28352, 2023. [58] L. Sadek, O. Sadek, H. T. Alaoui, M. S. Abdo, K. Shah, and T. Abdeljawad. Frac- S. Saber et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6152 24 of 24 tional order modeling of predicting covid-19 with isolation and vaccination strategies in morocco. CMES-Computer Modeling in Engineering & Sciences, 136:1931–1950, 2023. [59] O. Sadek, L. Sadek, S. Touhtouh, and A. Hajjaji. The mathematical fractional mod- eling of tio-2 nanopowder synthesis by sol–gel method at low temperature. Mathe- matical Modeling and Computing, 9(3):616–626, 2022. [60] H. P. Himsworth and S. K. Ker. The measurement of insulin sensitivity in normal and diabetic persons. Journal of Physiology, 1939. [61] V. W. Bolie. Coefficients of normal blood glucose regulation. Journal of Applied Physiology, 1961. [62] R. N. Bergman, Y. Z. Ider, C. R. Bowden, and C. Cobelli. Quantitative estimation of insulin sensitivity. American Journal of Physiology-Endocrinology and Metabolism, 1979. [63] R. N. Bergman, D. T. Finegood, and M. Ader. Assessment of insulin sensitivity in vivo. Endocrine Reviews, 1985. [64] M. Derouich and A. Boutayeb. The effect of physical exercise on the dynamics of glucose-insulin. Mathematical Biosciences, 2002. [65] G. Gaetano and O. Arino. Delay differential model of the glucose-insulin system. Mathematical Biosciences and Engineering, 2000. [66] M. Shabestari, R. Salimifar, and M. Rabiee. Chaotic behavior of glucose-insulin regulatory system using a predator-prey model. Chaos, Solitons & Fractals, 114:1– 10, 2018. [67] A. A. Elsadany. Complex dynamics in a fractional-order predator-prey model. Non- linear Dynamics, 67(4):2281–2289, 2012. [68] M. Caputo and M. Fabrizio. Application of new time and spatial fractional derivatives with exponential kernels. Progress in Fractional Differentiation and Applications, 2:1– 11, 2016. [69] J. Losada and J. J. Nieto. Properties of a new fractional derivative without singular kernel. Progress in Fractional Differentiation and Applications, 1:87–92, 2015. [70] A. Shaikh, A. Tassaddiq, K. S. Nisar, and D. Baleanu. Analysis of differential equa- tions involving caputo-fabrizio fractional operator and its applications to reaction- diffusion equations. Advances in Difference Equations, 2019(1):178, 2019. [71] S. A. Khan et al. Existence theory and numerical solutions to smoking model under caputo-fabrizio fractional derivative. Chaos: An Interdisciplinary Journal of Nonlin- ear Science, 29(1):013128, 2019.