Introduction

Chronic metabolic diseases like diabetes mellitus remain a major worldwide health concern. It is distinguished by chronic hyperglycemia brought on by insufficiency in either the production or the action of insulin, or both. These physiological abnormalities highlight the critical function of insulin as an anabolic hormone by interfering with the metabolism of proteins, fats, and sugar. Diabetics are at a higher risk of coronary artery disease and have a fourfold increased risk of stroke. Untreated diabetes can have dangerous consequences, such as vision problems that increase the risk of infection, blindness, and unconsciousness. However, certain people, especially children, may exhibit obvious symptoms such as impaired vision, increased appetite (polyphagia), excessive thirst (polydipsia), excessive urination (polyuria), and abrupt weight loss if they have a full insulin deficiency. The International Diabetes Federation (IDF) estimated that 415 million individuals aged 20 to 79 had diabetes mellitus in 2015. According to projections, population concern will have risen by an additional 200 million come the year 20401. Today, the problem is global, and according to the data obtained in 2021, about 537 million people, or 10.5% of the world’s population suffer from diabetes2. It is projected that these numbers could increase to 643 million by 2030 and increase to 783 million by 20453. These statistics bring out the importance of screening and assessment approaches in diabetes intervention. The etiopathogenesis of diabetes mellitus genetic factors comprises an extensive range including changes in genes concerning insulin synthesis, secretion, and degradation. Diabetes mellitus can be defined by multiple distortions in protein and lipid balance, as well as global glucose balance. The knowledge of this type of diabetes mellitus has several etiology, the detailed look at the different pathophysiologic processes involved, and developing treatment strategies to address those pathophysiologic processes has been the purpose of many researches as knowledge in this area expands4. There is a close relationship between COVID-19 and diabetes; COVID-19 worsens diabetes and makes its management more challenging. Some preliminary research shows that patients with diabetes who get infected with COVID-19 have worsened outcomes and have more extended symptoms, like fatigue and higher glucose levels. For example, one research found that 44% of diabetic patients tested for high blood sugar after being cured of COVID-19, and the long-term side effects of the virus are consistent with Chronic Fatigue Syndrome, such as fatigue and headache5. COVID-19 from a physiological perspective might affect the pancreas causing inflammation and damaging its function to provide adequate insulin leading to strong hyperglycemia complications6. Besides, COVID-19 has caused higher levels of stress and the disruption of normal diabetic care, worse glycemic control, and higher rates of incident diabetes. The Centers for Disease Control and Prevention has stated that patients with Type 1 and Type 2 diabetes are among those at a higher risk for the severity of COVID-19 due to impaired immune functions7. Thus, the pandemic does not only endanger the primary lives of diabetes patients but also comes with long-term implications for diabetes care and health results.

Atangana8 has introduced a completely innovative methodology for fractional calculus referred to as the fractal fractional derivative. The traditional methodology is augmented by this innovative technique about fractal fractions9,10. This development facilitates the simultaneous examination of fractional operators and fractal dimensions during the analysis of fractal fractional derivatives. Creating models that more precisely represent memory-effect systems is a huge advantage of this operator. Here some latest research related to fractional, the Mittag–Leffler function is used in a COVID-19 epidemic model11 to better understand the complicated dynamics and stability of disease dissemination. Its findings help in the development of more accurate estimates and control methods for epidemic control. The interplay between diabetes and COVID-19 using a fractal-fractional model12, offering insights into how co-infections impact disease dynamics. The study13 establishes conditions for the uniqueness and existence of positive solutions for a class of fractional integro-differential equations, highlighting how specific constraints ensure these solutions. Identifies conditions under which positive solutions exist or do not exist for a fractional coupled system involving the generalized p-Laplacian14, shedding light on the role of fractional operators and nonlinearity in solution behavior. Fractional capacities to bounded open Lipschitz sets and their complements15, providing insights into how these capacities behave under fractional order settings and their impact on boundary conditions in mathematical analysis. The study16 introduces the M-truncate derivative to define the coupled fractional nonlinear Hirota equation and employs fractional functional and simple equation methods to derive new periodic and solitary wave solutions, enhancing the understanding of such equations. In17, an innovative computational techniques for solving fractional coupled nonlinear Helmholtz equations, providing robust solutions and demonstrating their efficiency through numerical simulations. Traces the fractional temperature field18, revealing how fractional-order derivatives model the diffusion process, capturing anomalous heat transfer behaviors more accurately than classical models. 19provided an insight into its long-term fluctuations and variance scaling, which differ from standard Brownian motion due to memory and dependency structures. A highly efficient numerical scheme for fractional evolution equations of distinct types20, showcasing its accuracy and stability in various applications. The study21 proposes new optical soliton solutions for fractional perturbed Schrödinger equations, advancing the understanding of optical fiber communication systems. Ebola virus dynamics using a piecewise hybrid fractional operator with the Mittag–Leffler kernel22, enhancing the understanding of epidemic progression and control on diverse time scales. A mathematical model of malaria transmission in Yemen23, utilizing piecewise Caputo-Fabrizio derivatives to capture complex dynamics and predict cases from 2022 to 2024, aiding in disease control and prevention. A novel fractional mathematical model using a \(\Phi\)-piecewise hybrid fractional derivative approach to analyze Ebola virus disease (EVD) transmission dynamics24, providing valuable insights for developing effective control strategies. The study25 introduces a Hepatitis B virus (HBV) model incorporating asymptomatic carriers, utilizing piecewise fractional-order derivatives to capture complex transmission dynamics, and provides comprehensive theoretical and numerical analysis to enhance understanding of HBV spread. A dengue model with fractal-fractional derivatives analyzes dynamics26, stability, and computing the basic reproduction number to enhance understanding of dengue transmission and inform effective control strategies. In27, a fractional-order model for forest resource conservation, addressing chaos control in the context of toxin activity and human-caused fires. The study28 uses a fractional-order kinetic model to analyze oleic acid epoxidation, with a focus on Ulam–Hyers stability. Its insights help improve reaction stability and efficiency in chemical processes involving oleic acid epoxidation. In study29 applies computational techniques to monitor a fractional-order model of type-1 diabetes, aiding in the feedback design of artificial pancreas systems. In30, the glucose-insulin-glucagon regulatory system in type I diabetes is analyzed using the Mittag–Leffler kernel, providing insights into the intricate dynamics of blood sugar regulation.

Fractional calculus has a long history and is used widely in the simulation of real-world physical phenomena. Fractional calculus studies have attracted a lot of attention lately. In31,  the diabetes model study is given; we modified the existing model and developed new results. The fractal-fractional Mittag–Leffler function has become an essential tool in mathematical modeling, particularly in systems with complex, memory-dependent dynamics like biological processes. In diabetes research, where the body’s glucose-insulin interaction is both nonlinear and subject to fractional memory effects, this approach holds significant promise. The fractional order model offers a deeper understanding of glucose-insulin regulation. Capturing complex dynamics, empowers more accurate monitoring and management of diabetes. Embrace this cutting-edge approach to improve patient outcomes. Together, we can transform healthcare with precision and innovation. The fractal-fractional Mittag–Leffler function can capture intricate feedback mechanisms and time-varying interactions, offering a nuanced view of glucose fluctuations and insulin regulation. Modeling diabetes with this framework enables researchers to explore more accurate predictions and personalized treatment strategies, advancing the understanding and management of this chronic condition in ways that traditional models simply cannot achieve. The choice of the fractal fractional derivative and the Mittag–Leffler kernel is motivated by their ability to model memory effects and complex, non-local interactions in glucose-insulin dynamics. These operators capture the inherent irregularities and scale-dependent behavior of biological systems more effectively than classical models. The Mittag–Leffler kernel’s flexibility allows for precise control over the decay and persistence of past states, while the fractal fractional derivative accommodates self-similar patterns in the system. Together, they provide a more accurate and robust framework for modeling diabetes dynamics compared to traditional differential approaches. Existing integer-order models of diabetes often fail to capture the complex, non-linear interactions and memory effects inherent in glucose-insulin dynamics. These models assume instantaneous responses, neglecting the gradual influence of past states on the system. In contrast, the proposed fractional-order model accounts for these memory effects, offering a more realistic representation of glucose regulation and improved control over glycemic variability. The fractal fractional derivative operator is employed in this research to capture the complex, dynamic nature of glucose-insulin interactions more effectively than traditional approaches, ultimately enhancing the understanding and management of diabetes.

The following sections of the document are organized as described: “Basic concepts” presents several key definitions and notations related to fractional calculus. In “Formulation of model”, the formulation of a FF-order mathematical model using the Mittag–Leffler kernel is introduced. “Existence and uniqueness” analyzes the proposed model, showing the existence and uniqueness of the solution through fixed point theory. In “Ulam–Hyers stability”, we investigate the Ulam–Hyers Stability. “Chaos Control” focuses on controlling chaos. In “PID and controllability”, we evaluate PID and Controllability. “Numerical scheme” contains the numerical analysis of the model. In “Results of proposed scheme”, numerical simulations and graphical results for the examined system are provided, employing the Newton polynomial method. Finally, “Conclusion” wraps up our findings.

Basic concepts

Here, we provide some key definitions that may help improve system analysis.

Definition 2.1

32 The fractal-fractional derivative of order \(\xi\) for a fractal differentiable and continuous function \(\mathfrak {B}(t)\) on the interval (a, b) can be formulated in terms of the Riemann–Liouville derivative when \(0 \le \xi , \eta \le 1\).

$$\begin{aligned} _{0}^{FFP}\textrm{D}_{t}^{\xi ,\zeta } \mathfrak {B}(t) = \frac{1}{\Gamma (n-\xi )} \frac{d}{dt^\zeta }\int ^{t}_0 (t-\epsilon )^{n-\xi -1} \mathfrak {A}(\epsilon ) d\epsilon . \end{aligned}$$
(1)

Where \(n-1<\xi ,\zeta <n\in \mathfrak {N}\). And \(\frac{\textrm{D} \mathfrak {B}(\epsilon )}{\textrm{D}\epsilon ^\zeta }=\lim _{t\rightarrow \epsilon } \frac{\mathfrak {B}(t)-\mathfrak {B}(\epsilon )}{t^\xi -\epsilon ^\zeta }\)

$$\begin{aligned} _{0}^{FFE}\textrm{D}_{t}^{\xi ,\zeta } \mathfrak {B}(t) = \frac{\mathfrak {Y}(\xi )}{\Gamma (n-\xi )} \frac{d}{dt^\zeta }\int ^{t}_0 \exp \Big [-\frac{\xi }{1-\xi } (t-\epsilon )\Big ] \mathfrak {B}(\epsilon )d\epsilon \end{aligned}$$
(2)

where \(\xi >0, \zeta \le n \in \mathfrak {N}\), and \(\mathfrak {Y}(0)=1=\mathfrak {Y}(1)\).

$$\begin{aligned} _{0}^{FFM}\textrm{D}_{t}^{\xi ,\zeta } \mathfrak {B}(t) = \frac{\mathscr{A}\mathscr{B}(\xi )}{1-\xi } \frac{d}{dt^\zeta }\int ^{t}_0 \mathfrak {B}(\epsilon ) \mathscr {F}_{\xi } \Big [-\frac{\xi }{1-\xi } (t-\epsilon )^\xi \Big ] d\epsilon . \end{aligned}$$
(3)

Here, \(0 < \xi , \zeta \le 1\), \(\mathfrak {F}_{\xi }\) represents the Mittag-Leffler function, and \(\mathscr{A}\mathscr{B}(\xi ) = 1 - \xi + \frac{\xi }{\Gamma (\xi )}\) denotes a normalization function.

Definition 2.2

8 A continuous function \(\mathfrak {B}(t)\) on (a, b) has a fractional order \(\xi\) and a fractal dimension \(\zeta\) when \(0 \le \xi , \zeta \le 1\).

$$\begin{aligned} _0^{FFP} \textrm{I} ^{\xi ,\zeta } \mathfrak {B}(t)= \frac{1}{\Gamma (\xi )} \int ^t_0 (t-\epsilon )^{\xi -1} \epsilon ^{1-\zeta }\mathfrak {B}(\epsilon )d\epsilon . \end{aligned}$$
(4)
$$\begin{aligned} _0^{FFE} \textrm{I} ^{\xi ,\zeta } \mathfrak {B}(t)= \frac{\zeta (1-\xi ) t^{\zeta -1} \mathfrak {B}(t)}{\mathfrak {Y}(\xi )}+\frac{\xi \zeta }{\mathfrak {Y}(\xi )} \int ^t_0 \epsilon ^{\xi -1} \mathfrak {B}(\epsilon ) d\epsilon . \end{aligned}$$
(5)
$$\begin{aligned} _0^{FFM} \textrm{I} ^{\xi ,\zeta } \mathfrak {B}(t)= \frac{\zeta (1-\zeta ) t^{\xi -1} \mathfrak {B}(t)}{\mathscr{A}\mathscr{B}(\xi )}+\frac{\xi \zeta }{\mathscr{A}\mathscr{B}(\xi )\Gamma (\xi )} \int ^t_0 (t-\epsilon )^{\xi -1} \epsilon ^{\zeta -1} \mathfrak {B}(\epsilon ) d\epsilon . \end{aligned}$$
(6)

Formulation of model

The dynamics of \(\beta\)-cell \((\beta ),\) glucose (G), insulin (I), and glucagon (J) are the main emphasis of the diabetes model incorporating parameters \(\beta\),I,G, and J are study in31. The model captures the interplay between \(\beta\)-cells, insulin, glucose, and glucagon in regulating blood sugar levels. Parameters like \(\sigma\), d and e determine how effectively glucose is regulated. The balance between insulin activity \(\sigma\), glucagon production p, and therapy b influences whether the system remains in a healthy state or progresses toward diabetes. Understanding these parameters helps in developing targeted treatments for diabetes management. The parameters values are shown in Table 1.

The Mittag–Leffler kernel is a generalization of the exponential function, commonly used in fractional calculus to describe processes with memory and hereditary properties. It provides a flexible framework for modeling systems where the effects of past states gradually diminish over time, unlike the constant decay in classical models. Fractal fractional operators combine fractional calculus with fractal geometry, extending the concept of differentiation and integration to account for dynamics on irregular or self-similar geometries. These operators are particularly useful in capturing the complexity of biological systems, such as glucose-insulin regulation, where memory effects and heterogeneous interactions play a significant role.

This diabetes model provides a comprehensive framework for understanding and managing glucose-insulin dynamics through advanced fractional techniques.

$$\begin{aligned} \left\{ \begin{array}{l} _0^{FFM}D_t^{\xi ,\varsigma }\beta \left( t \right) = \Lambda G + lJ - \left( {\sigma I + \rho + \nu } \right) \beta \\ _0^{FFM}D_t^{\xi ,\varsigma }I\left( t \right) = \frac{{\sigma \beta I}}{{e + G}} - \left( {d + \nu } \right) I\\ _0^{FFM}D_t^{\xi ,\varsigma }G\left( t \right) = bJ - \left( {\rho + \nu } \right) G\\ _0^{FFM}D_t^{\xi ,\varsigma }J\left( t \right) = p + \rho G - \left( {l + \nu + b} \right) J\end{array} \right. \end{aligned}$$
(7)

in terms of the initial conditions

$$\begin{aligned} \beta (0)\ge 0 \, \ \ \ I(0)\ge 0 \, \ \ \ G(0)\ge 0 \, \ \ \ J(0)\ge 0. \end{aligned}$$
Table 1 Model parameters and their biological relevance.
Full size table

Positively invariant region

We demonstrate the enclosed or closed set in this section.

$$\begin{aligned} \bigwedge =\{ \beta (t), I (t), G (t), J (t) \in \mathbb {R}_{+}^4\}. \end{aligned}$$
(8)

The region derived from the closed set, which serves as the system (7) feasibility region and guarantees positive invariance, is displayed in this section.

Lemma 1

The region \(\bigwedge \in \mathbb {R}_{+}^4\)

$$\begin{aligned} \bigwedge =\{( \beta (t), I (t), G (t), J (t))\in \mathbb {R}_{+}^4\}. \end{aligned}$$
(9)

It draws all solutions from the system (7) and maintains positive invariance, given non-negative initial conditions within the proposed system in \(\mathbb {R}_{+}^4\).

Proof

Here are the outcomes of our depiction of the system’s (7) positive solutions:

$$\begin{aligned} \begin{aligned}&{_0^{FFM}D_t^{\xi ,\varsigma }\beta \left( t \right) \left| {_{\beta = 0}} \right. = \Lambda G + lJ \ge 0,}\\ &{_0^{FFM}D_t^{\xi ,\varsigma }I\left( t \right) \left| {_{I = 0}} \right. = 0,}\\ &{_0^{FFM}D_t^{\xi ,\varsigma }G\left( t \right) \left| {_{G = 0}} \right. = bJ \ge 0,}\\ &{_0^{FFM}D_t^{\xi ,\varsigma }J\left( t \right) \left| {_{J = 0}} \right. = p + \rho G \ge 0} \end{aligned} \end{aligned}$$
(10)

The vector field stays inside the region \(\mathbb {R}_{+}^4\) along each hyperplane containing the non-negative orthant for \(t \ge 0\), according to the system (10). Therefore, we can say that the area \(\bigwedge\),

$$\begin{aligned} \bigwedge =\{( \beta , I (t), G (t), J (t))\in \mathbb {R}_{+}^4\}, \end{aligned}$$

constitutes a region that maintains positive invariance. \(\square\)

Equilibrium point analysis

The equilibrium point for the proposed model is as follows:

$$\begin{aligned} \begin{aligned} {E^ \star }&= \left( {{\beta _ \star },{I_ \star },{G_ \star },{J_ \star }} \right) \\&\left( {\frac{{p(\Lambda b + l\nu + l\rho )}}{{(\nu + \rho )(b\nu + l\nu + l\rho + \nu \rho + {\nu ^2})}},0,\frac{{bp}}{{b\nu + l\nu + l\rho + \nu \rho + {\nu ^2}}},\frac{{\nu p + p\rho }}{{b\nu + l\nu + l\rho + \nu \rho + {\nu ^2}}}} \right) \end{aligned} \end{aligned}$$
(11)

Existence and uniqueness

In this section, we elucidate the existence and uniqueness of the solution for the proposed model through the application of the fixed-point methodology.

$$\begin{aligned} \left\{ {\begin{array}{l}{_0^{FFM}D_t^{\xi ,\varsigma }\beta \left( t \right) = \Lambda G + lJ - \left( {\sigma I + \rho + \nu } \right) \beta = \beta \left( {t,\hat{\beta }\left( t \right) } \right) }\\ {_0^{FFM}D_t^{\xi ,\varsigma }I\left( t \right) = \frac{{\sigma \beta I}}{{e + G}} - \left( {d + \nu } \right) I = I\left( {t,{\hat{I}}\left( t \right) } \right) }\\ {_0^{FFM}D_t^{\xi ,\varsigma }G\left( t \right) = bJ - \left( {\rho + \nu } \right) G = G\left( {t,{\hat{G}}\left( t \right) } \right) }\\ {_0^{FFM}D_t^{\xi ,\varsigma }J\left( t \right) = p + \rho G - \left( {l + \nu + b} \right) J = J\left( {t,{\hat{J}}\left( t \right) } \right) }\end{array}} \right. \end{aligned}$$
(12)

Expression reformulated using Fractal-Fractional integral and Mittag–Leffler kernel.

$$\begin{aligned} \begin{aligned}&\beta (t)=\beta (0)+\frac{\varsigma (1-\xi ) t^{\varsigma -1}}{{A B}(\xi )} \hat{\beta }\left( t, \beta (t)\right) +\frac{\xi \varsigma }{{A B}(\xi ) \Gamma (\xi )} \int _0^t(t-\epsilon )^{\xi -1} \epsilon ^{\varsigma -1} \hat{\beta }\left( \epsilon , \beta (\epsilon )\right) d \epsilon =\mathscr {A}_1+\mathscr {A}_2, \\&I (t)= I (0)+\frac{\varsigma (1-\xi ) t^{\varsigma -1}}{{A B}(\xi )} \hat{ I }\left( t, I (t)\right) +\frac{\xi \varsigma }{{A B}(\xi ) \Gamma (\xi )} \int _0^t(t-\epsilon )^{\xi -1} \epsilon ^{\varsigma -1} \hat{ I }\left( \epsilon , I (\epsilon )\right) d \epsilon =\mathscr {B}_1+\mathscr {B}_2, \\&G (t)= G (0)+\frac{\varsigma (1-\xi ) t^{\varsigma -1}}{{A B}(\xi )} \hat{ G }\left( t, G (t)\right) +\frac{\xi \varsigma }{{A B}(\xi ) \Gamma (\xi )} \int _0^t(t-\epsilon )^{\xi -1} \epsilon ^{\varsigma -1} \hat{ G }\left( \epsilon , G (\epsilon )\right) d \epsilon =\mathscr {C}_1+\mathscr {C}_2, \\&J (t)= J (0)+\frac{\varsigma (1-\xi ) t^{\varsigma -1}}{{A B}(\xi )} \hat{ J }\left( t, J (t)\right) +\frac{\xi \varsigma }{{A B}(\xi ) \Gamma (\xi )} \int _0^t(t-\epsilon )^{\xi -1} \epsilon ^{\varsigma -1} \hat{ J }\left( \epsilon , J (\epsilon )\right) d \xi =\mathscr {D}_1+\mathscr {D}_2. \end{aligned} \end{aligned}$$
(13)

We show that \(\zeta (\mathscr {A}_1, \mathscr {B}_1, \mathscr {C}_1, \mathscr {D}_1)\), the principal elements of the governing equations (13), operate as contraction maps. Moreover, \(\varpi (\mathscr {A}_2, \mathscr {B}_2, \mathscr {C}_2, \mathscr {D}_2)\) functions as continuous compact integral component parts. The fixed point theorem of Krasnoselski is used for these validations.

$$\begin{aligned} \begin{aligned}&{\mathscr {A}_1} = \beta (0) + \frac{{\varsigma (1 - \xi ){t^{\varsigma - 1}}}}{{AB(\xi )}}\hat{\beta }\left( {t,\beta (t)} \right) ,{\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\mathscr {A}_2} = \frac{{\xi \varsigma }}{{AB(\xi )\Gamma (\xi )}}\int _0^t {{{(t - \epsilon )}^{\xi - 1}}} {\epsilon ^{\varsigma - 1}}\hat{\beta }\left( {\epsilon ,\beta (\epsilon )} \right) d\epsilon ,\\ &{\mathscr {B}_1} = I(0) + \frac{{\varsigma (1 - \xi ){t^{\varsigma - 1}}}}{{AB(\xi )}}{\hat{I}}\left( {t,I(t)} \right) ,{\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\mathscr {B}_2} = \frac{{\xi \varsigma }}{{AB(\xi )\Gamma (\xi )}}\int _0^t {{{(t - \epsilon )}^{\xi - 1}}} {\epsilon ^{\varsigma - 1}}{\hat{I}}\left( {\xi ,I(\xi )} \right) d\epsilon ,\\ &{\mathscr {C}_1} = G(0) + \frac{{\varsigma (1 - \xi ){t^{\varsigma - 1}}}}{{AB(\xi )}}{\hat{G}}\left( {t,G(t)} \right) ,{\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\mathscr {C}_2} = \frac{{\xi \varsigma }}{{AB(\xi )\Gamma (\xi )}}\int _0^t {{{(t - \epsilon )}^{\xi - 1}}} {\epsilon ^{\varsigma - 1}}{\hat{G}}\left( {\epsilon ,G(\epsilon )} \right) d\epsilon ,\\ &{\mathscr {D}_1} = J(0) + \frac{{\varsigma (1 - \xi ){t^{\varsigma - 1}}}}{{AB(\xi )}}{\hat{J}}\left( {t,J(t)} \right) ,{\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\hspace{0.55542pt}} {\mathscr {D}_2} = \frac{{\xi \varsigma }}{{AB(\xi )\Gamma (\xi )}}\int _0^t {{{(t - \xi )}^{\xi - 1}}} {\epsilon ^{\varsigma - 1}}{\hat{J}}\left( {\epsilon ,J(\epsilon )} \right) d\epsilon ,\end{aligned} \end{aligned}$$
(14)

Theorem 4.1

The transformation \(\zeta (\mathscr {A}_1, \mathscr {B}_1, \mathscr {C}_1, \mathscr {D}_1):[0, \mathbb {T}]\rightarrow R^4\) is nonlinear. When \(\beth _{\rho ^{1}}, \beth _{\rho ^{2}}, \beth _{\rho ^{3}}, \beth _{\rho ^{4}}>0\), as stated in equation (13), guarantees a Lipschitz contractive condition.

Proof

A completely normed space contains the definition of the operator \(\zeta (\mathscr {A}_1, \mathscr {B}_1, \mathscr {C}_1, \mathscr {D}_1):[0, \mathbb {T}]\rightarrow R^4\). The norm here denotes

$$\begin{aligned} \begin{aligned} |\beta , I , G , J |=&\max _{t \in [0, \mathbb {T}]}|| \beta (t), I (t), G (t), J (t)||\\&\text {where}~~ \beta , I , G , J \in [0, \mathbb {T}]. \end{aligned} \end{aligned}$$
(15)

(i) First, we shall show that \(\zeta (\mathscr {A}_1, \mathscr {B}_1, \mathscr {C}_1, \mathscr {D}_1)\) functions as a contraction map. Taking into account \(\beta (t)\) and \(\hat{\beta }(t)\), we find that

$$\begin{aligned} \begin{aligned}&\left\| {{\rho ^1}(\beta ,I,G,J)(t) - {\rho ^1}(\hat{\beta },I,G,J)(t)} \right\| \\ &= \left\| {\left( {\Lambda G + lJ - \left( {\sigma I + \rho + \nu } \right) \beta } \right) - \left( {\Lambda G + lJ - \left( {\sigma I + \rho + \nu } \right) \hat{\beta }} \right) } \right\| \\ &= \left\| {\left( {\sigma \left\| I \right\| + \rho + \nu } \right) \left( {\beta - \hat{\beta }} \right) } \right\| \\ &\le \left\| {\left( {\sigma \left\| I \right\| + \rho + \nu } \right) } \right\| \left\| {\left( {\beta - \hat{\beta }} \right) } \right\| \le {\beth _{{\rho ^1}}}\left\| {\left( {\beta - \hat{\beta }} \right) } \right\| , \end{aligned} \end{aligned}$$
(16)

\(\left\| {{\rho ^1}(\beta ,I,G,J)(t) - {\rho ^1}(\hat{\beta },I,G,J)(t)} \right\| \le {\beth _{{\rho ^1}}}\left\| {\left( {\beta - \hat{\beta }} \right) } \right\|\)   where  \({\beth _{{\rho ^1}}} = \left\| {\left( {\sigma \left\| I \right\| + \rho + \nu } \right) } \right\| .\)

Similarly for (I, G, J),

\(\left\| {\rho ^{2}\left( {\beta ,I,G,J} \right) - \rho ^{2}\left( {\beta ,{\hat{I}},G,J} \right) } \right\| \le {\beth _{\rho ^2} }\left\| {I - {\hat{I}}} \right\|\)   where  \({\beth _{\rho ^2} } = \left\| {\left( {\frac{{\sigma \left\| \beta \right\| }}{{e + \left\| G \right\| }} - \left( {d + \nu } \right) } \right) } \right\| ,\)

\(\left\| {\rho ^{3}\left( {\beta ,I,G,J} \right) - \rho ^{3}\left( {\beta ,I,{\hat{G}},J} \right) } \right\| \le {\beth _{\rho ^3} }\left\| {G - {\hat{G}}} \right\|\)   where  \({\beth _{\rho ^3} } = \left\| {\left( {\rho + \nu } \right) } \right\| ,\)

\(\left\| {\rho ^{4}\left( {\beta ,I,G,J} \right) - \rho ^{4}\left( {\beta ,I,G,{\hat{J}}} \right) } \right\| \le {\beth _{\rho ^4} }\left\| {J - {\hat{J}}} \right\|\)   where  \({\beth _{\rho ^4} } = \left\| {\left( {l + \nu + b} \right) } \right\| ,\)

This indicates that in the case of the operator \(\zeta ( \beta , I , G , J )\), we observe

$$\begin{aligned} \begin{aligned}&||\zeta ( \beta I , G , J )(t)-\zeta (\hat{\beta },\hat{I},\hat{G},\hat{J})|| \\&\quad =\Psi \max _{t \in [0, \mathbb {T}]}||(\beta , I , G , J )(t)-(\hat{\beta },\hat{I},\hat{G},\hat{J})(t)||\\&\quad \le \Psi ||(\beta , I , G , J )(t)-(\hat{\beta },\hat{I},\hat{G},\hat{J})(t)|| \\&\quad \le \Psi \beth , \end{aligned} \end{aligned}$$
(17)

where \(\Psi = \frac{\varsigma (1-\xi ) t^{\varsigma -1}}{{A B}(\xi )}\). With \(\beth =\max [\beth _{\rho ^{1}}, \beth _{\rho ^{2}}, \beth _{\rho ^{3}}, \beth _{\rho ^{4}}]<1\) being a Lipschitz constant, it follows that \(\varsigma (\rho ^1, \rho ^2, \rho ^3, \rho ^4)\) represents a non-expansive operator.

To prove \(\varpi (\mathscr {A}_2, \mathscr {B}_2, \mathscr {C}_2, \mathscr {D}_2)\) is Continuously compact:

We will demonstrate the continuity and compactness of \(\varpi (\mathscr {A}_2, \mathscr {B}_2, \mathscr {C}_2, \mathscr {D}_2)\). The absolute value of all positively bounded continuous operators \(\rho ^1, \rho ^2, \rho ^3, \rho ^4\) specified in equations (13), indicated by non-zero positive constants \(\mathscr {Q}_{\rho ^{1}}, \mathscr {Q}_{\rho ^{2}}, \mathscr {Q}_{\rho ^{3}}, \mathscr {Q}_{\rho ^{4}}\), \(\mathscr {Y}_{\rho ^{1}}, \mathscr {Y}_{\rho ^{2}}, \mathscr {Y}_{\rho ^{3}}, \mathscr {Y}_{\rho ^{4}}\), satisfying certain boundedness conditions, confirm the compactness of the operator \(\varpi (\rho ^{1}, \rho ^{2}, \rho ^{3}, \rho ^{4})\).

$$\begin{aligned} \begin{aligned}&|\rho ^{1}(t,\beta )| \le \mathscr {Q}_{\rho ^{1}}||\beta ||+\mathscr {Y}_{\rho ^{1}}, \\&|\rho ^{2}(t, I )| \le \mathscr {Q}_{\rho ^{2}}|| I ||+\mathscr {Y}_{\rho ^{2}}, \\&|\rho ^{3}(t, G )| \le \mathscr {Q}_{\rho ^{3}}|| G ||+\mathscr {Y}_{\rho ^{3}}, \\&|\rho ^{4}(t, J )| \le \mathscr {Q}_{\rho ^{4}}|| J ||+\mathscr {Y}_{\rho ^{4}}. \end{aligned} \end{aligned}$$
(18)

If we consider \(\mathfrak {P}\) to be a closed subset of \(\mathfrak {B}\), then

$$\begin{aligned} \mathfrak {P}=\{(\rho ^1, \rho ^2, \rho ^3, \rho ^4) \in \mathfrak {B} /||\rho ^1, \rho ^2, \rho ^3, \rho ^4|| \le \varphi , \varphi >0\}. \end{aligned}$$
(19)

When \((\rho ^1, \rho ^2, \rho ^3, \rho ^4) \in \mathfrak {P}\), we get

$$\begin{aligned} \begin{aligned} \left\| {\mathscr {A}_2}\left( t,\beta \right) \right\|&=\max _{t \in [0, \mathbb {T}]}\left| \frac{\xi \varsigma }{{AB}(\xi ) \Gamma (\xi )} \int _0^t(t-\epsilon )^{\xi -1} \epsilon ^{\varsigma -1} \rho ^1\left( \epsilon , \beta (\epsilon )\right) d \epsilon \right| \\&\le \Psi _1 \int _0^{\eta }(\eta -\epsilon )^{\xi -1} \epsilon ^{\varsigma -1}\left| \rho ^1\left( \epsilon , \beta (\epsilon )\right) \right| d \epsilon \\&\le \Psi _1 \mathscr {Q}_{\rho ^{1}} \varphi +\mathscr {Y}_{\rho ^{1}}. \end{aligned} \end{aligned}$$
(20)

where \(\Psi _1 ={\frac{{\xi \varsigma \Gamma \left( \varsigma \right) }}{{{AB}\left( \xi \right) \Gamma \left( {\xi + \varsigma } \right) }}}\). In this manner, we have

$$\begin{aligned} \begin{aligned}&\left\| {{\mathscr {B}_2}\left( {t,I} \right) } \right\| \le \Psi _1 \mathscr {Q}_{\rho ^{2}}\varphi + \mathscr {Y}_{\rho ^{2}},\\ &\left\| {{\mathscr {C}_2}\left( {t,G} \right) } \right\| \le \Psi _1 \mathscr {Q}_{\rho ^{3}}\varphi + \mathscr {Y}_{\rho ^{3}},\\ &\left\| {{\mathscr {D}_2}\left( {t,J} \right) } \right\| \le \Psi _1 \mathscr {Q}_{\rho ^{4}}\varphi + \mathscr {Y}_{\rho ^{4}}.\end{aligned} \end{aligned}$$

For the other compartments in the suggested model, we use a similar method to determine that the maximal norm of \(\Vert \varpi (\mathscr {A}_2, \mathscr {B}_2, \mathscr {C}_2, \mathscr {D}_2) \Vert\) is:

$$\begin{aligned} \begin{aligned}&\Vert \varpi (\mathscr {A}_2, \mathscr {B}_2, \mathscr {C}_2, \mathscr {D}_2)\Vert \\&\le \Psi _1[(\mathscr {Q}_{\rho ^{1}}+\mathscr {Q}_{\rho ^{2}}+\mathscr {Q}_{\rho ^{3}}+\mathscr {Q}_{\rho ^{4}}) \varphi + \mathscr {Y}_{\rho ^{1}}+ \mathscr {Y}_{\rho ^{2}}+ \mathscr {Y}_{\rho ^{3}}+ \mathscr {Y}_{\rho ^{4}}]=\Im . \end{aligned} \end{aligned}$$
(21)

We prove that \(\Vert \mathscr {H}(\beta , I , G , J ) \Vert \le \Im\), where \(\Im\) is a positive constant. This suggests that \(\mathscr {H}\) is a uniformly bounded operator.

Next, we will demonstrate the equi-continuity of \(\mathscr {H}\) for \(t_a < t_b \in [0, \mathbb {T}]\). To achieve this, for \(t_a < t_b \in [0, \mathbb {T}]\), we have

$$\begin{aligned} \left\{ \begin{array}{l} \left| \mathscr {A}_{2}\left( t_b, \beta \right) -\mathscr {A}_{2}\left( t_a, \beta \right) \right| \\ \quad =\frac{\xi \varsigma }{{A B}(\xi ) \Gamma (\xi )}\left| \int _0^{t_b}(t-\epsilon )^{\xi -1} \epsilon ^{\varsigma -1} \mathscr {A}_{2}\left( \epsilon , \beta (\epsilon )\right) d \epsilon -\int _0^{t_a}(t-\epsilon )^{\xi -1} \epsilon ^{\varsigma -1} \mathscr {A}_{2} \left[ \epsilon , \beta (\epsilon )\right] d \epsilon \right| \\ \quad \le \frac{\xi \varsigma }{{AB}(\xi ) \Gamma (\xi )}\left[ \int _0^{t_b}(t-\epsilon )^{\xi -1} \epsilon ^{\varsigma -1}-\int _0^{t_a}(t-\epsilon )^{\xi -1} \epsilon ^{\varsigma -1}\right] \left( \mathscr {Q}_{\rho ^{1}} \varphi +\mathscr {Y}_{\rho ^{1}}\right) \\ \quad \le [\Psi _1 \mathscr {Q}_{\rho ^{1}} \varphi +\mathscr {Y}_{\rho ^{1}}]\left( t_b-t_a\right) ^{\xi +\varsigma -1}. \end{array}\right. \end{aligned}$$
(22)

Likewise, in a similar manner

$$\begin{aligned} \begin{aligned}&\left| \mathscr {B}_{2}\left( t_b, I \right) -\mathscr {B}_{2}\left( t_a, I \right) \right| \le [\Psi _1 \mathscr {Q}_{\rho ^{2}} \varphi +\mathscr {Y}_{\rho ^{2}}]\left( t_b-t_a\right) ^{\xi +\varsigma -1}, \\&\left| \mathscr {C}_{2}\left( t_b, G \right) -\mathscr {C}_{2}\left( t_a, G \right) \right| \le [\Psi _1 \mathscr {Q}_{\rho ^{3}} \varphi +\mathscr {Y}_{\rho ^{3}}]\left( t_b-t_a\right) ^{\xi +\varsigma -1},\\&\left| \mathscr {D}_{2}\left( t_b, J \right) -\mathscr {D}_{2}\left( t_a, J \right) \right| \le [\Psi _1 \mathscr {Q}_{\rho ^{4}} \varphi +\mathscr {Y}_{\rho ^{4}}]\left( t_b-t_a\right) ^{\xi +\varsigma -1}. \end{aligned} \end{aligned}$$
(23)

Given that \(t_b \rightarrow t_a\) is unrelated to \((\beta , I , G , J )\), this infers that

$$\begin{aligned} \left\| \mathscr {H}(\beta , I , G , J )\left( t_2\right) -\mathscr {H}(\beta , I , G , J )\left( t_1\right) \right\| \rightarrow 0. \end{aligned}$$
(24)

Consequently, \(\mathscr {H}(\mathscr {A}_{2}, \mathscr {B}_{2}, \mathscr {C}_{2}, \mathscr {D}_{2})\) denotes an operator that is both totally continuous and equi-continuous. \(\mathscr {H}(\mathscr {A}_{2}, \mathscr {B}_{2}, \mathscr {C}_{2}, \mathscr {D}_{2})\) is comparatively compact according to Arzela’s theorem. According to Krasnoselski’s theorem, the operators \(\zeta\) and \(\varpi\) guarantee the existence of a unique solution due to their contraction and continuity qualities.\(\square\)

Theorem 4.2

The model (7) possesses a unique solution if

$$\begin{aligned} \left[ \Psi _1 + \Psi \right] \beth \le 1, \end{aligned}$$

where \(\beth =\max \{\beth _{\rho ^{1}}, \beth _{\rho ^{2}}, \beth _{\rho ^{3}}, \beth _{\rho ^{4}}\}\).

Proof

Define an operator \(\Pi =(\Pi _1, \Pi _2, \Pi _3, \Pi _4): \mathfrak {B} \rightarrow \mathfrak {B}\) using Eq. (13) in the following manner:

$$\begin{aligned} \left\{ \begin{array}{l} \Pi _1\left( \beta \right) (t)=\beta (0)+\Psi \rho ^{1}\left( t, \beta (t)\right) +\frac{\xi \varsigma }{{A B}(\xi ) \Gamma (\xi )} \int _0^t(t-\epsilon )^{\xi -1} \epsilon ^{\varsigma -1} \rho ^{1}\left( \epsilon , \beta (\epsilon )\right) d \epsilon , \\ \Pi _2\left( I \right) (t)= I (0)+\Psi \rho ^{2}\left( t, I (t)\right) +\frac{\xi \varsigma }{{A B}(\xi ) \Gamma (\xi )} \int _0^t(t-\epsilon )^{\xi -1} \epsilon ^{\varsigma -1} \rho ^{2}\left( \epsilon , I (\epsilon )\right) d \epsilon , \\ \Pi _3\left( G \right) (t)= G (0)+\Psi \rho ^{3}\left( t, G (t)\right) +\frac{\xi \varsigma }{{A B}(\xi ) \Gamma (\xi )} \int _0^t(t-\epsilon )^{\xi -1} \epsilon ^{\varsigma -1} \rho ^{3}\left( \epsilon , G (\epsilon )\right) d \epsilon , \\ \Pi _4\left( J \right) (t)= J (0)+\Psi \rho ^{4}\left( t, J (t)\right) +\frac{\xi \varsigma }{{A B}(\xi ) \Gamma (\xi )} \int _0^t(t-\epsilon )^{\xi -1} \epsilon ^{\varsigma -1} \rho ^{4}\left( \epsilon , J (\epsilon )\right) d \epsilon . \\ \end{array}\right. \end{aligned}$$

Considering \((\beta , I , G , J ), (\hat{\beta }, \hat{ I }, \hat{ G }, \hat{ J }) \in \mathfrak {B}\) and employing equation (20), we obtain

$$\begin{aligned} \begin{aligned} \left\| \Pi _1\left( \beta \right) (t)-\Pi _1\left( \hat{\beta }\right) (t)\right\|&=\Psi \left\| \rho ^{1}\left( t, \beta (t)\right) -\rho ^{1}\left( t, \hat{\beta }(t)\right) \right\| \\&+\frac{\xi \varsigma }{{AB}(\xi ) \Gamma (\xi )} \int _0^t\left\| \rho ^{1}\left( \epsilon , \beta (\epsilon )\right) -\rho ^{1}\left( \epsilon , \hat{\beta }(\epsilon )\right) \right\| (t-\epsilon )^{\xi -1} \epsilon ^{\varsigma -1} \rho ^{1}\left( \epsilon , \beta (\epsilon )\right) d \epsilon \\&\le \Psi \beth _{\rho ^{1}}\left\| \beta -\hat{\beta }\right\| +\Psi _1 \beth _{\rho ^{1}}\left\| \beta -\hat{\beta }\right\| \\&\le \left[ \Psi +\Psi _1\right] \beth _{\rho ^{1}}\left\| \beta -\hat{\beta }\right\| , \end{aligned} \end{aligned}$$

\(\left\| \beta -\hat{\beta }\right\| \rightarrow 0\) when \(\beta \rightarrow \hat{\beta }\). Hence

$$\begin{aligned} \left\| \Pi _1\left( \beta \right) (t)-\Pi _1\left( \hat{\beta }\right) (t)\right\| \le \left[ \Psi + \Psi _1\right] \beth _{\rho ^{1}} \le 1, \end{aligned}$$

with

$$\begin{aligned} \left\| \Pi _1\left( \beta \right) (t)-\Pi _1\left( \hat{\beta }\right) (t)\right\| \left[ 1-\left( \Psi +\Psi _1\right) \beth _{\rho ^{1}}\right] \le 0. \end{aligned}$$

Continuing this method, we determine

$$\begin{aligned} \begin{aligned}&\left\| \Pi \left( \beta , I , G , J \right) (t)-\Pi \left( \hat{\beta }, \hat{ I }, \hat{ G }, \hat{ J }\right) (t)\right\| \\&\le \left( \Psi + \Psi _1\right) \beth \left\| \left( \beta , I , G , J \right) -\left( \hat{\beta }, \hat{ I }, \hat{ G }, \hat{ J }\right) \right\| . \end{aligned} \end{aligned}$$

\(\square\)

Ulam–Hyers stability

$$\begin{aligned} \begin{aligned}&\left| {\beta \left( t \right) + \left( {\frac{{\varsigma \left( {1 - \xi } \right) {t^{\varsigma - 1}}}}{{{AB}\left( \xi \right) }}} \right) {o _1}\left( {t,\beta \left( t \right) } \right) + (\frac{{\xi \varsigma }}{{{AB}\left( \xi \right) }}\int \limits _0^t {\left( {t - \epsilon } \right) {\epsilon ^{\xi - 1}}){o _1}\left( {t,\beta \left( t \right) } \right) d\epsilon } } \right| \le \P _1,\\ &\left| {I\left( t \right) + \left( {\frac{{\varsigma \left( {1 - \xi } \right) {t^{\varsigma - 1}}}}{{{AB}\left( \xi \right) }}} \right) {o _2}\left( {t,I \left( t \right) } \right) + (\frac{{\xi \varsigma }}{{{AB}\left( \xi \right) }}\int \limits _0^t {\left( {t - \epsilon } \right) {\epsilon ^{\xi - 1}}){o _2}\left( {t,I \left( t \right) } \right) d\epsilon } } \right| \le \P _2,\\ &\left| {G\left( t \right) + \left( {\frac{{\varsigma \left( {1 - \xi } \right) {t^{\varsigma - 1}}}}{{{AB}\left( \xi \right) }}} \right) {o _3}\left( {t,G \left( t \right) } \right) + (\frac{{\xi \varsigma }}{{{AB}\left( \xi \right) }}\int \limits _0^t {\left( {t - \epsilon } \right) {\epsilon ^{\xi - 1}}){o _3}\left( {t,G \left( t \right) } \right) d\epsilon } } \right| \le \P _3,\\ &\left| {J\left( t \right) + \left( {\frac{{\varsigma \left( {1 - \xi } \right) {t^{\varsigma - 1}}}}{{{AB}\left( \xi \right) }}} \right) {o _4}\left( {t,J \left( t \right) } \right) + (\frac{{\xi \varsigma }}{{{AB}\left( \xi \right) }}\int \limits _0^t {\left( {t - \epsilon } \right) {\epsilon ^{\xi - 1}}){o _3}\left( {t,J \left( t \right) } \right) d\epsilon } } \right| \le \P _4.\end{aligned} \end{aligned}$$
(25)

there exist \({\hat{\beta }}, {{\hat{I}}}, {{\hat{G}}}, {{\hat{J}}}\) are satisfying

$$\begin{aligned} \begin{aligned}&\hat{\beta }\left( t \right) = \left( {\frac{{\varsigma \left( {1 - \xi } \right) {t^{\varsigma - 1}}}}{{{A}{B}\left( \xi \right) }}} \right) {o _1}\left( {t,\hat{\beta }\left( t \right) } \right) + (\frac{{\xi \varsigma }}{{{A}{B}\left( \xi \right) }}\int \limits _0^t {\left( {t - \epsilon } \right) {\epsilon ^{\xi - 1}}){o _1}\left( {t,\hat{\beta }\left( t \right) } \right) d\epsilon },\\ &{\hat{I}}\left( t \right) = \left( {\frac{{\varsigma \left( {1 - \xi } \right) {t^{\varsigma - 1}}}}{{{A}{B}\left( \xi \right) }}} \right) {o _2}\left( {t,{\hat{I}} \left( t \right) } \right) + (\frac{{\xi \varsigma }}{{{A}{B}\left( \xi \right) }}\int \limits _0^t {\left( {t - \epsilon } \right) {\epsilon ^{\xi - 1}}){o _2}\left( {t,{\hat{I}} \left( t \right) } \right) d\epsilon }, \\ &{\hat{G}}\left( t \right) = \left( {\frac{{\varsigma \left( {1 - \xi } \right) {t^{\varsigma - 1}}}}{{{A}{B}\left( \xi \right) }}} \right) {o _3}\left( {t,{\hat{G}} \left( t \right) } \right) + (\frac{{\xi \varsigma }}{{{A}{B}\left( \xi \right) }}\int \limits _0^t {\left( {t - \epsilon } \right) {\epsilon ^{\xi - 1}}){o _3}\left( {t,{\hat{G}} \left( t \right) } \right) d\epsilon }, \\ &{\hat{J}}\left( t \right) = \left( {\frac{{\varsigma \left( {1 - \xi } \right) {t^{\varsigma - 1}}}}{{{A}{B}\left( \xi \right) }}} \right) {o _4}\left( {t,{\hat{J}}\left( t \right) } \right) + (\frac{{\xi \varsigma }}{{{A}{B}\left( \xi \right) }}\int \limits _0^t {\left( {t - \epsilon } \right) {\epsilon ^{\xi - 1}}){o _4}\left( {t,{\hat{J}} \left( t \right) } \right) d\epsilon }\end{aligned} \end{aligned}$$
(26)

such that

$$\begin{aligned} \beta \left( t \right) - {\hat{G}} \le o_1 \Upsilon _1, \end{aligned}$$

Similarly for other

$$\begin{aligned} \begin{aligned}&I\left( t \right) - {\hat{I}} \le o_2 \Upsilon _2,\\ &G\left( t \right) - {\hat{G}} \le o_3 \Upsilon _3,\\ &J\left( t \right) - {\hat{J}} \le o_4 \Upsilon _4. \end{aligned} \end{aligned}$$
(27)

Theorem 5.1

The functions \(\beta \left( t \right) ,I \left( t \right) ,G \left( t \right) , J \left( t \right)\) and \(\hat{\beta }\left( t \right) ,{\hat{I}}\left( t \right) ,{\hat{G}}\left( t \right) , {\hat{J}}\left( t \right)\) are consistent under certain conditions

$$\begin{aligned} \left\| {\beta \left( t \right) } \right\| \le {\wp _1},\left\| {I \left( t \right) } \right\| \le {\wp _2},\left\| {G \left( t \right) } \right\| \le {\wp _3},\left\| {J \left( t \right) } \right\| \le {\wp _4} \end{aligned}$$
(28)

Ulam–Hyers stability for the fractional model (7) is achieved if a specific condition is satisfied.

Proof

A solution has been determined for \(\beta (t)\), \(I(t)\), \(G(t)\), and \(J(t)\). Let \(\hat{\beta }(t)\), \(\hat{I}(t)\), \(\hat{G}(t)\), and \(\hat{J}(t)\) be approximate solutions that satisfy the specified conditions for the system.

$$\begin{aligned} \begin{aligned} \left\| {\beta \left( t \right) - \hat{\beta }\left( t \right) } \right\|&\le \frac{{\varsigma \left( {1 - \xi } \right) {t^{\varsigma - 1}}}}{{{AB}\left( \xi \right) }}{\wp _1}\left( {t,\hat{\beta }\left( t \right) } \right) + \frac{{\xi \varsigma }}{{{AB}\left( \xi \right) }}\int \limits _0^t {{{\left( {t - \epsilon } \right) }^{\xi - 1}}{\wp _1}\left( {t,\hat{\beta }\left( t \right) } \right) d\epsilon } \\ &\le \frac{{\varsigma \left( {1 - \xi } \right) {t^{\varsigma - 1}} + \xi \varsigma }}{{{AB}\left( \xi \right) }}{\wp _1}\left\| {\beta \left( t \right) - \hat{\beta }\left( t \right) } \right\| \end{aligned} \end{aligned}$$
(29)

consider \({\wp _i} = o_i\) and \(\Upsilon _i = \frac{{\varsigma \left( {1 - \xi } \right) {t^{\varsigma - 1}} + \xi \varsigma }}{{{AB}\left( \xi \right) }}\) for \(i = 1,2,3,4\)

$$\begin{aligned} \left\| {G\left( t \right) - {\bar{G}}\left( t \right) } \right\| \le \varpi _1 \Omega _1 \end{aligned}$$
(30)

similarly,

$$\begin{aligned} \begin{aligned} \left\| {I \left( t \right) - {\hat{I}} \left( t \right) } \right\|&\le \frac{{\varsigma \left( {1 - \xi } \right) {t^{\varsigma - 1}}}}{{{AB}\left( \xi \right) }}{\wp _2}\left( {t,{\hat{I}} \left( t \right) } \right) + \frac{{\xi \varsigma }}{{{AB}\left( \xi \right) }}\int \limits _0^t {{{\left( {t - \epsilon } \right) }^{\xi - 1}}{\wp _2}\left( {t,{\hat{I}} \left( t \right) } \right) d\epsilon } \\ &\le \frac{{\varsigma \left( {1 - \xi } \right) {t^{\varsigma - 1}} + \xi \varsigma }}{{{AB}\left( \xi \right) }}{\wp _2}\left\| {I \left( t \right) - {\hat{I}} \left( t \right) } \right\| \end{aligned} \end{aligned}$$
(31)
$$\begin{aligned} \begin{aligned} \left\| {G \left( t \right) - {\hat{G}} \left( t \right) } \right\|&\le \frac{{\varsigma \left( {1 - \xi } \right) {t^{\varsigma - 1}}}}{{{AB}\left( \xi \right) }}{\wp _3}\left( {t,{\hat{G}} \left( t \right) } \right) + \frac{{\xi \varsigma }}{{{AB}\left( \xi \right) }}\int \limits _0^t {{{\left( {t - \epsilon } \right) }^{\xi - 1}}{\wp _3}\left( {t,{\hat{G}} \left( t \right) } \right) d\epsilon } \\ &\le \frac{{\varsigma \left( {1 - \xi } \right) {t^{\varsigma - 1}} + \xi \varsigma }}{{{AB}\left( \xi \right) }}{\wp _3}\left\| {G \left( t \right) - {\hat{G}} \left( t \right) } \right\| \end{aligned} \end{aligned}$$
(32)

and

$$\begin{aligned} \begin{aligned} \left\| {J \left( t \right) - {\hat{J}} \left( t \right) } \right\|&\le \frac{{\varsigma \left( {1 - \xi } \right) {t^{\varsigma - 1}}}}{{{AB}\left( \xi \right) }}{\wp _4}\left( {t,{\hat{J}} \left( t \right) } \right) + \frac{{\xi \varsigma }}{{{AB}\left( \xi \right) }}\int \limits _0^t {{{\left( {t - \epsilon } \right) }^{\xi - 1}}{\wp _4}\left( {t,{\hat{J}} \left( t \right) } \right) d\epsilon } \\ &\le \frac{{\varsigma \left( {1 - \xi } \right) {t^{\varsigma - 1}} + \xi \varsigma }}{{{AB}\left( \xi \right) }}{\wp _4}\left\| {J \left( t \right) - {\hat{J}} \left( t \right) } \right\| \end{aligned} \end{aligned}$$
(33)

These contrasts demonstrate the system Ulam-Hyers stability. \(\square\)

Thus, the connection of Ulam–Hyers stability to the proposed glucose-insulin model of fractional order is justified by the fact that this condition prevents major changes in the function behavior due to small perturbations in the control space. This property ensures that the model’s forecasts do not greatly fluctuate, or deviate in unpredicted ways, given uncertainties or changes in patient information. Thus, proving the Ulam-Hyers stability, the model can be confidently used for monitoring changes in glucose levels in real conditions and improving its use in clinical practice for diabetic patients.

Chaos control

In this section, we use the linear feedback control method to stabilize system (7) to equilibrium points. The fractional-order system in its controlled form (7) will be analyzed as follows:

$$\begin{aligned} \left\{ {\begin{array}{l}{_0^{FFM}D_t^{\xi ,\varsigma }\beta \left( t \right) = \Lambda G + lJ - \left( {\sigma I + \rho + \nu } \right) \beta + {\mu _1}\left( {{\beta _t} - \beta } \right) }\\ {_0^{FFM}D_t^{\xi ,\varsigma }I\left( t \right) = \frac{{\sigma \beta I}}{{e + G}} - \left( {d + \nu } \right) I + {\mu _2}\left( {{I_t} - I} \right) }\\ {_0^{FFM}D_t^{\xi ,\varsigma }G\left( t \right) = bJ - \left( {\rho + \nu } \right) G + {\mu _3}\left( {{G_t} - G} \right) }\\ {_0^{FFM}D_t^{\xi ,\varsigma }J\left( t \right) = p + \rho G - \left( {l + \nu + b} \right) J + {\mu _4}\left( {{J_t} - J} \right) }\end{array}} \right. \end{aligned}$$
(34)

where \(E^\star\) (7) is the system equilibrium point and the control variables are \(\mu _1\), \(\mu _2\), \(\mu _3\), and \(\mu _4\). For the system (34), the Jacobian matrix at \(E^\star\) is built as

$$\begin{aligned} J(E^\star ) = \begin{bmatrix} - \nu - \rho - \mu _1 & -\frac{p(\Lambda b + l \nu + l \rho )}{(\nu + \rho )(b \nu + l \nu + l \rho + \nu \rho + \nu ^2)} \sigma & \Lambda & l \\ 0 & \frac{\frac{p(\Lambda b + l \nu + l \rho )}{(\nu + \rho )(b \nu + l \nu + l \rho + \nu \rho + \nu ^2)} \sigma }{\frac{b p}{b \nu + l \nu + l \rho + \nu \rho + \nu ^2} + e} - \nu - \mu _2 - d & 0 & 0 \\ 0 & 0 & - \nu - \rho - \mu _3 & b \\ 0 & 0 & \rho & - b - l - \nu - \mu _4 \\ \end{bmatrix} \end{aligned}$$

\(\mu _1=1,\mu _2=2,\mu _3=3,\mu _4=4\) is the selection. The typical roots of the polynomial with respect to the equilibrium point \(E^\star\) are as follows.

$$\begin{aligned} C\left( \lambda \right) = {\lambda ^4} +85.3236{\lambda ^3} + 48.4057{\lambda ^2} + 11.6145\lambda + 53.6499. \end{aligned}$$
(35)

Moreover, \(\lambda _1= -1.6100\), \(\lambda _2= -4.2460\), \(\lambda _3= -3.5440\), and \(\lambda _4= -2.2145\) are computed as the characteristic roots of Eq. (35). Given that each of the eigenvalues is a negative real number, \(E^\star\) is an asymptotically stable equilibrium point. Surface phase plots in 3D at different fractional order values of the GIG model derived in Figs. 1 and 2 at different initial conditions show the results in feasible control. We can easily observe the effect of the relation of beta cells and glucagon effect on glucose-insulin compartments to manage the glucose level in the human body.

Chaotic behavior in diabetes dynamics suggests unpredictable fluctuations in glucose and insulin levels, complicating disease management. This highlights the need for robust, adaptive therapeutic strategies like real-time glucose monitoring and dynamic insulin delivery systems. Understanding chaos improves prediction models, aiding in early detection of destabilizing patterns. These insights pave the way for personalized and precision-based interventions to optimize patient outcomes.

Fig. 1
figure 1

GIG Surface phase plot at different fractional order values in feasible control.

Full size image
Fig. 2
figure 2

GIG Surface phase plot at different fractional order and initial values in feasible control.

Full size image

PID and controllability

The diabetic mathematical model is provided as

$$\begin{aligned} \left\{ {\begin{array}{l}{\dot{\beta }\left( t \right) = \Lambda G + lJ - \left( {\sigma I + \rho + \nu } \right) \beta }\\ {\dot{I}\left( t \right) = \frac{{\sigma \beta I}}{{e + G}} - \left( {d + \nu } \right) I}\\ {\dot{G}\left( t \right) = bJ - \left( {\rho + \nu } \right) G}\\ {\dot{J}\left( t \right) = p + \rho G - \left( {l + \nu + b} \right) J}\end{array}} \right. \end{aligned}$$
(36)

In this study, we use the Mittag–Leffler kernel to build a fractal-fractional model for diabetes. The following is a description of this model:

$$\begin{aligned} \left\{ \begin{array}{l} _0^{FFM}D_t^{\xi ,\varsigma }\beta \left( t \right) = \Lambda G + lJ - \left( {\sigma I + \rho + \nu } \right) \beta \\ _0^{FFM}D_t^{\xi ,\varsigma }I\left( t \right) = \frac{{\sigma \beta I}}{{e + G}} - \left( {d + \nu } \right) I\\ _0^{FFM}D_t^{\xi ,\varsigma }G\left( t \right) = bJ - \left( {\rho + \nu } \right) G\\ _0^{FFM}D_t^{\xi ,\varsigma }J\left( t \right) = p + \rho G - \left( {l + \nu + b} \right) J\end{array} \right. \end{aligned}$$
(37)

where \(0 < \xi \le 1\).

Fractional-order PID controller33 is as follows:

$$\begin{aligned} \begin{aligned}&{v_1}\left( t \right) = {q_p}{c_1}\left( t \right) + {q_i}_0\gimel _t^{ - \xi }{c_1}\left( t \right) + {q_d}^{FFM}D_{0,t}^{\xi ,\varsigma }{c_1}\left( t \right) ,\\ &{v_2}\left( t \right) = {q_p}{c_2}\left( t \right) + {q_i}_0\gimel _t^{ - \xi }{c_2}\left( t \right) + {q_d}^{FFM}D_{0,t}^{\xi ,\varsigma }{c_2}\left( t \right) ,\\ &{v_3}\left( t \right) = {q_p}{c_3}\left( t \right) + {q_i}_0\gimel _t^{ - \xi }{c_3}\left( t \right) + {q_d}^{FFM}D_{0,t}^{\xi ,\varsigma }{c_3}\left( t \right) ,\\ &{v_4}\left( t \right) = {q_p}{c_4}\left( t \right) + {q_i}_0\gimel _t^{ - \xi }{c_4}\left( t \right) + {q_d}^{FFM}D_{0,t}^{\xi ,\varsigma }{c_4}\left( t \right) ,\end{aligned} \end{aligned}$$
(38)

where \(c_1\left( t \right) = \left( {\beta \left( t \right) - {\beta ^*}} \right) ,\,\,c_2\left( t \right) = \left( {I\left( t \right) - {I^*}} \right)\), \(c_3\left( t \right) = \left( {G \left( t \right) - {G ^*}} \right)\)  and  \(c_4\left( t \right) = \left( {J \left( t \right) - J^*} \right)\) and \(q_{p},q_{i}\) and \(q_d\) are the control gains.

The system controller (38) includes proportional, integral, and derivative components governed by gains \(q_p\), \(q_i\), and \(q_d\). It simplifies to a fractional PD controller when \(q_i = 0\), a fractional PI controller when \(q_d = 0\), and a conventional PID controller for \(\varepsilon = 1\).

Applying controller (38) to model (37) yields the following controlled model:

$$\begin{aligned} \begin{aligned}&{_0^{FFM}D_t^{\xi ,\varsigma }\beta \left( t \right) = \Lambda G\left( {t - \tau } \right) + lJ\left( {t - \tau } \right) - \left( {\sigma I\left( {t - \tau } \right) + \rho + \nu } \right) \beta \left( {t - \tau } \right) ,}\\ &{_0^{FFM}D_t^{\xi ,\varsigma }I\left( t \right) = \frac{{\sigma \beta \left( {t - \tau } \right) I\left( {t - \tau } \right) }}{{e + G\left( {t - \tau } \right) }} - \left( {d + \nu } \right) I\left( {t - \tau } \right) ,}\\ &{_0^{FFM}D_t^{\xi ,\varsigma }G\left( t \right) = bJ\left( {t - \tau } \right) - \left( {\rho + \nu } \right) G\left( {t - \tau } \right) ,}\\ &{_0^{FFM}D_t^{\xi ,\varsigma }J\left( t \right) = p + \rho G\left( {t - \tau } \right) - \left( {l + \nu + b} \right) J\left( {t - \tau } \right) .}\end{aligned} \end{aligned}$$
(39)

Declare the recently added variables.

$$\begin{aligned} \begin{aligned}&{\daleth _1}\left( t \right) { = _0}\gimel _t^{ - \xi }{c_1}\left( t \right) ,\\ &{\daleth _2}\left( t \right) { = _0}\gimel _t^{ - \xi }{c_2}\left( t \right) ,\\ &{\daleth _3}\left( t \right) { = _0}\gimel _t^{ - \xi }{c_3}\left( t \right) ,\\ &{\daleth _4}\left( t \right) { = _0}\gimel _t^{ - \xi }{c_4}\left( t \right) ,\end{aligned} \\ \begin{aligned}&_0^{FFM}D_t^{\xi ,\varsigma }\gimel _1\left( t \right) = c_1\left( t \right) ,\\&_0^{FFM}D_t^{\xi ,\varsigma }\daleth _2\left( t \right) = c_2\left( t \right) ,\\&_0^{FFM}D_t^{\xi ,\varsigma }\daleth _3\left( t \right) = c_3\left( t \right) ,\\&_0^{FFM}D_t^{\xi ,\varsigma }\daleth _4\left( t \right) = c_4\left( t \right) .\end{aligned} \end{aligned}$$

As thus, the fractional-order controlled model (39) can be expressed as

$$\begin{aligned} \begin{aligned}&{_0^{FFM}D_t^{\xi ,\varsigma }\beta \left( t \right) = \left[ {\Lambda G\left( {t - \tau } \right) + lJ\left( {t - \tau } \right) - \left( {\sigma I\left( {t - \tau } \right) + \rho + \nu } \right) \beta \left( {t - \tau } \right) + {q_p}\left( {\beta \left( t \right) - {\beta ^ * }} \right) + {q_i}{\daleth _1}} \right] \frac{1}{{1 - {q_d}}},}\\ &{_0^{FFM}D_t^{\xi ,\varsigma }I\left( t \right) = \left[ {\frac{{\sigma \beta \left( {t - \tau } \right) I\left( {t - \tau } \right) }}{{e + G\left( {t - \tau } \right) }} - \left( {d + \nu } \right) I\left( {t - \tau } \right) + {q_p}\left( {I\left( t \right) - {I^ * }} \right) + {q_i}{\daleth _2}} \right] \frac{1}{{1 - {q_d}}},}\\ &{_0^{FFM}D_t^{\xi ,\varsigma }G\left( t \right) = \left[ {bJ\left( {t - \tau } \right) - \left( {\rho + \nu } \right) G\left( {t - \tau } \right) + {q_p}\left( {G\left( t \right) - {G^ * }} \right) + {q_i}{\daleth _3}} \right] \frac{1}{{1 - {q_d}}},}\\ &{_0^{FFM}D_t^{\xi ,\varsigma }J\left( t \right) = \left[ {p + \rho G\left( {t - \tau } \right) - \left( {l + \nu + b} \right) J\left( {t - \tau } \right) + {q_p}\left( {J\left( t \right) - {J^ * }} \right) + {q_i}{\daleth _4}} \right] \frac{1}{{1 - {q_d}}}.}\end{aligned} \end{aligned}$$
(40)
$$\begin{aligned} \begin{aligned}&_0^{FFM}D_t^{\xi ,\varsigma }{\daleth _1}\left( t \right) = \left( {\beta \left( t \right) - {\beta ^*}} \right) ,\\ &_0^{FFM}D_t^{\xi ,\varsigma }{\daleth _2}\left( t \right) = \left( {{I}\left( t \right) - I^ * } \right) ,\\ &_0^{FFM}D_t^{\xi ,\varsigma }{\daleth _3}\left( t \right) = \left( {{G}\left( t \right) - G^ * } \right) ,\\ &_0^{FFM}D_t^{\xi ,\varsigma }{\daleth _4}\left( t \right) = \left( {J\left( t \right) - {J^ * }} \right) .\end{aligned} \end{aligned}$$

Analyzing the controlled model (39) is challenging due to integral variable terms in the controller (38). To simplify, variable transformation is applied, yielding the corresponding-order model (40) for easier analysis.

The controlled model (40) has a unique positive equilibrium, where the first components \(\beta ^\bullet\), \(I^\bullet\), \(G^\bullet\), and \(J^\bullet\) match the equilibrium of the uncontrolled model (38). This shows the fractional-order PID controller (39) effectively preserves the original model unique equilibrium.

Numerical scheme

In this part, we numerically solve the system (7) using the Newton polynomial.

$$\begin{aligned} \left\{ \begin{array}{l} _0^{FFM}D_t^{\xi ,\varsigma }\beta \left( t \right) = \Lambda G + lJ - \left( {\sigma I + \rho + \nu } \right) \beta \\ _0^{FFM}D_t^{\xi ,\varsigma }I\left( t \right) = \frac{{\sigma \beta I}}{{e + G}} - \left( {d + \nu } \right) I\\ _0^{FFM}D_t^{\xi ,\varsigma }G\left( t \right) = bJ - \left( {\rho + \nu } \right) G\\ _0^{FFM}D_t^{\xi ,\varsigma }J\left( t \right) = p + \rho G - \left( {l + \nu + b} \right) J\end{array} \right. \end{aligned}$$
(41)

To make things simple, we could describe the system mentioned above as follows:

$$\begin{aligned} \begin{aligned}& ^{FFM}\textrm{D}_{0,t}^{\xi ,\varsigma }\beta (t)= \wp _1[t,\beta ,I,G,J],\\ & ^{FFM}\textrm{D}_{0,t}^{\xi ,\varsigma }I(t)= \wp _2[t,\beta ,I,G,J],\\ & ^{FFM}\textrm{D}_{0,t}^{\xi ,\varsigma }G(t)= \wp _3[t,\beta ,I,G,J],\\ & ^{FFM}\textrm{D}_{0,t}^{\xi ,\varsigma }J(t)= \wp _4[t,\beta ,I,G,J]. \end{aligned} \end{aligned}$$
(42)

Using the Mittag–Leffler kernel and the fractal-fractional integral, we get the following:

$$\begin{aligned} \beta (t_{b} +1)&= \beta (0)+\frac{1-\xi }{{AB}(\xi )} t_{b}^{1-\varsigma } \wp _1\big [t_{b},\beta (t_{b}),I(t_{b}),G(t_{b}),J(t_{b}))\big ] \nonumber \\&+ \frac{\xi }{{AB}(\xi )\Gamma (\xi )}\displaystyle \sum _{e=2}^{b} \int _{t_e}^{t_{e+1}} \wp _1(t,\beta ,I,G,J)\epsilon ^{1-\varsigma }(t_{{b}+1}-\epsilon )^{\xi -1}d\epsilon \end{aligned}$$
(43)
$$\begin{aligned} I(t_{b} +1)&= I(0)+\frac{1-\xi }{{AB}(\xi )} t_{b}^{1-\varsigma } \wp _2\big [t_{b},\beta (t_{b}),I(t_{b}),G(t_{b}),J(t_{b}))\big ] \nonumber \\&+ \frac{\xi }{{AB}(\xi )\Gamma (\xi )}\displaystyle \sum _{e=2}^{b} \int _{t_e}^{t_{e+1}} \wp _2(t,\beta ,I,G,J)\epsilon ^{1-\varsigma }(t_{{b}+1}-\epsilon )^{\xi -1}d\epsilon \end{aligned}$$
(44)
$$\begin{aligned} G(t_{b} +1)&= G(0)+\frac{1-\xi }{{AB}(\xi )} t_{b}^{1-\varsigma } \wp _3\big [t_{b},\beta (t_{b}),I(t_{b}),G(t_{b}),J(t_{b}))\big ] \nonumber \\&+ \frac{\xi }{{AB}(\xi )\Gamma (\xi )}\displaystyle \sum _{e=2}^{b} \int _{t_e}^{t_{e+1}} \wp _3(t,\beta ,I,G,J)\epsilon ^{1-\varsigma }(t_{{b}+1}-\epsilon )^{\xi -1}d\epsilon \end{aligned}$$
(45)
$$\begin{aligned} J(t_{b} +1)&= J(0)+\frac{1-\xi }{{AB}(\xi )} t_{b}^{1-\varsigma } \wp _4\big [t_{b},\beta (t_{b}),I(t_{b}),G(t_{b}),J(t_{b}))\big ] \nonumber \\&+ \frac{\xi }{{AB}(\xi )\Gamma (\xi )}\displaystyle \sum _{e=2}^{b} \int _{t_e}^{t_{e+1}} \wp _4(t,\beta ,I,G,J)\epsilon ^{1-\varsigma }(t_{{b}+1}-\epsilon )^{\xi -1}d\epsilon \end{aligned}$$
(46)

Here, the Newton polynomial is reviewed:

$$\begin{aligned} \mathfrak {B}(t,\beta ,I,G,J)&\simeq ~ \mathfrak {B}(t_{b-2},\beta _{b-2},I{_{b-2}},G{_{b-2}},J_{b-2})+\frac{1}{\Delta t}\Big [\mathfrak {B}(t_{b-1},\beta _{b-1},I{_{b-1}},G{_{b-1}},J_{b-1})\nonumber \\&-\mathfrak {B}(t_{b-2},\beta _{b-2},I{_{b-2}},G{_{b-2}},J_{b-2})\Big ] \times (\epsilon -t_{{b}-2})+\frac{1}{2{\Delta t}^2}\Big [\mathfrak {B}(t_{b},\beta _{b},I{_{b}},G{_{b}},J_{b})\nonumber \\&-2\mathfrak {B}(t_{b-1},\beta _{b-1},I{_{b-1}},G{_{b-1}},J_{b-1}) -\mathfrak {B}(t_{b-2},\beta _{b-2},I{_{b-2}},G{_{b-2}},J_{b-2})\Big ]\nonumber \\&\times (\epsilon -t_{{b}-2})(\epsilon -t_{{b}-1}) \end{aligned}$$
(47)

The Newton polynomial (47) can be used to solve equations (43)–(46). This provides us with

$$\begin{aligned} \begin{aligned} \beta _{(b +1)}&= \beta (0)+ \frac{1-\xi }{{AB}(\xi )} t_{b}^{1-\varsigma } \wp _1\big (t_{b},\beta (t_{b}),I(t_{b}), G(t_{b}),J(t_{b})\big )\\ &+\frac{\xi }{{AB}(\xi )\Gamma (\xi )}\sum _{e=2}^{b}\wp _1\big (t_{e-2},\beta ^{e-2},I^{e-2}, G^{e-2},J^{e-2}\big )t_{e-2}^{1-\varsigma }\times \int _{t_e}^{t_{e+1}}(t_{{b}+1}-\epsilon )^{\xi -1}d\epsilon \\ &+\frac{\xi }{{AB}(\xi )\Gamma (\xi )}\sum _{e=2}^{b} \frac{1}{\Delta t} \Big [t_{e-1}^{1-\varsigma }~~\wp _1\big (t_{e-1},\beta ^{e-1},I^{e-1},G^{e-1},J^{e-1}\big )\\ &-t_{e-2}^{1-\varsigma }~~\wp _1\big (t_{e-2},\beta ^{e-2},I^{e-2},G^{e-2},J^{e-2})\Big ] \times \int _{t_e}^{t_{e+1}} (\epsilon -t_{e-2})(t_{{b}+1}-\epsilon )^{\xi -1}d\epsilon \\ &+\frac{\xi }{{AB}(\xi )\Gamma (\xi )}\sum _{e=2}^{b} \frac{1}{2\Delta t^2} \Big [t_e^{1-\varsigma }\wp _1\big (t_{e},\beta ^{e},I^{e},G^{e},J^{e}\big ) \\ &-2t_{e-1}^{1-\varsigma }\wp _1\big (t_{e-1},\beta ^{e-1},I^{e-1},G^{e-1},J^{e-1}\big )\\ &+t_{e-2}^{1-\varsigma }\wp _1\big (t_{e-2},\beta ^{e-2},I^{e-2},G^{e-2},J^{e-2}\big )\Big ]\times \int _{t_e}^{t_{e+1}} (\epsilon -t_{e-2})(\epsilon -t_{e-1})(t_{{b}+1}-\epsilon )^{\xi -1}d\epsilon \end{aligned} \end{aligned}$$
(48)
$$\begin{aligned} \begin{aligned} I_{(b +1)}&= I(0)+ \frac{1-\xi }{{AB}(\xi )} t_{b}^{1-\varsigma } \wp _2\big (t_{b},\beta (t_{b}),I(t_{b}), G(t_{b}),J(t_{b})\big )\\ &+\frac{\xi }{{AB}(\xi )\Gamma (\xi )}\sum _{e=2}^{b}\wp _2\big (t_{e-2}, \beta ^{e-2},I^{e-2},G^{e-2},J^{e-2}\big )t_{e-2}^{1-\varsigma }\times \int _{t_e}^{t_{e+1}}(t_{{b}+1}-\epsilon )^{\xi -1}d\epsilon \\ &+\frac{\xi }{{AB}(\xi )\Gamma (\xi )}\sum _{e=2}^{b} \frac{1}{\Delta t} \Big [t_{e-1}^{1-\varsigma }~~\wp _2\big (t_{e-1}, \beta ^{e-1},I^{e-1},G^{e-1},J^{e-1}\big )\\ &-t_{e-2}^{1-\varsigma }~~\wp _2\big (t_{e-2},\beta ^{e-2},I^{e-2}, G^{e-2},J^{e-2})\Big ] \times \int _{t_e}^{t_{e+1}} (\epsilon -t_{e-2})(t_{{b}+1}-\epsilon )^{\xi -1}d\epsilon \\ &+\frac{\xi }{{AB}(\xi )\Gamma (\xi )}\sum _{e=2}^{b} \frac{1}{2\Delta t^2} \Big [t_e^{1-\varsigma }\wp _2\big (t_{e}, \beta ^{e},I^{e},G^{e},J^{e}\big ) \\ &-2t_{e-1}^{1-\varsigma }\wp _2\big (t_{e-1},\beta ^{e-1},I^{e-1},G^{e-1},J^{e-1}\big )\\ &+t_{e-2}^{1-\varsigma }\wp _2\big (t_{e-2},\beta ^{e-2},I^{e-2},G^{e-2},J^{e-2}\big )\Big ]\times \int _{t_e}^{t_{e+1}} (\epsilon -t_{e-2})(\epsilon -t_{e-1})(t_{{b}+1}-\epsilon )^{\xi -1}d\epsilon \end{aligned} \end{aligned}$$
(49)
$$\begin{aligned} \begin{aligned} G_{(b +1)}&= G(0)+ \frac{1-\xi }{{AB}(\xi )} t_{b}^{1-\varsigma } \wp _3\big (t_{b},\beta (t_{b}),I(t_{b}), G(t_{b}),J(t_{b})\big )\\ &+\frac{\xi }{{AB}(\xi )\Gamma (\xi )}\sum _{e=2}^{b}\wp _3\big (t_{e-2}, \beta ^{e-2},I^{e-2},G^{e-2},J^{e-2}\big )t_{e-2}^{1-\varsigma }\times \int _{t_e}^{t_{e+1}}(t_{{b}+1}-\epsilon )^{\xi -1}d\epsilon \\ &+\frac{\xi }{{AB}(\xi )\Gamma (\xi )}\sum _{e=2}^{b} \frac{1}{\Delta t} \Big [t_{e-1}^{1-\varsigma }~~\wp _3 \big (t_{e-1},\beta ^{e-1},I^{e-1},G^{e-1},J^{e-1}\big )\\ &-t_{e-2}^{1-\varsigma }~~\wp _3\big (t_{e-2}, \beta ^{e-2},I^{e-2},G^{e-2},J^{e-2})\Big ] \times \int _{t_e}^{t_{e+1}} (\epsilon -t_{e-2})(t_{{b}+1}-\epsilon )^{\xi -1}d\epsilon \\ &+\frac{\xi }{{AB}(\xi )\Gamma (\xi )}\sum _{e=2}^{b} \frac{1}{2\Delta t^2} \Big [t_e^{1-\varsigma }\wp _3 \big (t_{e},\beta ^{e},I^{e},G^{e},J^{e}\big ) \\ &-2t_{e-1}^{1-\varsigma }\wp _3\big (t_{e-1},\beta ^{e-1},I^{e-1},G^{e-1},J^{e-1}\big )\\ &+t_{e-2}^{1-\varsigma }\wp _3\big (t_{e-2},\beta ^{e-2},I^{e-2},G^{e-2},J^{e-2}\big )\Big ]\times \int _{t_e}^{t_{e+1}} (\epsilon -t_{e-2})(\epsilon -t_{e-1})(t_{{b}+1}-\epsilon )^{\xi -1}d\epsilon \end{aligned} \end{aligned}$$
(50)
$$\begin{aligned} \begin{aligned} J_{(b +1)}&= J(0)+ \frac{1-\xi }{{AB}(\xi )} t_{b}^{1-\varsigma } \wp _4\big (t_{b},\beta (t_{b}), I(t_{b}), G(t_{b}),J(t_{b})\big )\\ &+\frac{\xi }{{AB}(\xi )\Gamma (\xi )}\sum _{e=2}^{b}\wp _4\big (t_{e-2}, \beta ^{e-2},I^{e-2},G^{e-2},J^{e-2}\big )t_{e-2}^{1-\varsigma }\times \int _{t_e}^{t_{e+1}}(t_{{b}+1}-\epsilon )^{\xi -1}d\epsilon \\ &+\frac{\xi }{{AB}(\xi )\Gamma (\xi )}\sum _{e=2}^{b} \frac{1}{\Delta t} \Big [t_{e-1}^{1-\varsigma }~~\wp _4 \big (t_{e-1},\beta ^{e-1},I^{e-1},G^{e-1},J^{e-1}\big )\\ &-t_{e-2}^{1-\varsigma }~~\wp _4\big (t_{e-2}, \beta ^{e-2},I^{e-2},G^{e-2},J^{e-2})\Big ] \times \int _{t_e}^{t_{e+1}} (\epsilon -t_{e-2})(t_{{b}+1}-\epsilon )^{\xi -1}d\epsilon \\ &+\frac{\xi }{{AB}(\xi )\Gamma (\xi )}\sum _{e=2}^{b} \frac{1}{2\Delta t^2} \Big [t_e^{1-\varsigma }\wp _4 \big (t_{e},\beta ^{e},I^{e},G^{e},J^{e}\big ) \\ &-2t_{e-1}^{1-\varsigma }\wp _4\big (t_{e-1},\beta ^{e-1},I^{e-1},G^{e-1},J^{e-1}\big )\\ &+t_{e-2}^{1-\varsigma }\wp _4\big (t_{e-2},\beta ^{e-2},I^{e-2},G^{e-2},J^{e-2}\big )\Big ]\times \int _{t_e}^{t_{e+1}} (\epsilon -t_{e-2})(\epsilon -t_{e-1})(t_{{b}+1}-\epsilon )^{\xi -1}d\epsilon \end{aligned} \end{aligned}$$
(51)

we get

$$\begin{aligned} \begin{aligned} \begin{aligned}\beta (b + 1)&= \beta (0) + \frac{{1 - \xi }}{{AB(\xi )}}t_b^{1 - \varsigma }{\wp _1}[{t_b},\beta ({t_b}),I({t_b}),G({t_b}),J({t_b})]\\ &+ \frac{{\xi {{(\Delta t)}^\xi }}}{{AB(\xi )\Gamma (\xi + 1)}}\sum \limits _{e = 2}^b {{\wp _1}} [{t_{e - 2}},{\beta ^{e - 2}},{I^{e - 2}},{G^{e - 2}},{J^{e - 2}}]t_{e - 2}^{1 - \varsigma } \times \Xi \\ &+ \frac{{\xi {{(\Delta t)}^\xi }}}{{AB(\xi )\Gamma (\xi + 2)}}\sum \limits _{e = 2}^b \{ t_{e - 1}^{1 - \varsigma }{\wp _1}[{t_{e - 1}},{\beta ^{e - 1}},{I^{e - 1}},{G^{e - 1}},{J^{e - 1}}] - t_{e - 2}^{1 - \varsigma }{\wp _1}[{t_{e - 2}},{\beta ^{e - 2}},{I^{e - 2}},{G^{e - 2}},{J^{e - 2}}]\} \times {\Xi _1}\\ &+ \frac{{\xi {{(\Delta t)}^\xi }}}{{2AB(\xi )\Gamma (\xi + 3)}}\sum \limits _{e = 2}^b \{ t_e^{1 - \varsigma }{\wp _1}[{t_e},{\beta ^e},{I^e},{G^e},{J^e}] - 2t_{e - 1}^{1 - \varsigma }{\wp _1}[{t_{e - 1}},{\beta ^{e - 1}},{I^{e - 1}},{G^{e - 1}},{J^{e - 1}}]\\ &+ t_{e - 2}^{1 - \varsigma }{\wp _1}[{t_{e - 2}},{\beta ^{e - 2}},{I^{e - 2}},{G^{e - 2}},{J^{e - 2}}]\} \times {\Xi _2},\end{aligned} \end{aligned} \end{aligned}$$
(52)
$$\begin{aligned} \begin{aligned} I(b + 1)&= I(0) + \frac{{1 - \xi }}{{AB(\xi )}}t_b^{1 - \varsigma }{\wp _2}[{t_b},\beta ({t_b}),I({t_b}),G({t_b}),J({t_b})]\\ &+ \frac{{\xi {{(\Delta t)}^\xi }}}{{AB(\xi )\Gamma (\xi + 1)}}\sum \limits _{e = 2}^b {{\wp _2}} [{t_{e - 2}},{\beta ^{e - 2}},{I^{e - 2}},{G^{e - 2}},{J^{e - 2}}]t_{e - 2}^{1 - \varsigma } \times \Xi \\ &+ \frac{{\xi {{(\Delta t)}^\xi }}}{{AB(\xi )\Gamma (\xi + 2)}}\sum \limits _{e = 2}^k \{ t_{e - 1}^{1 - \varsigma }{\wp _2}[{t_{e - 1}},{\beta ^{e - 1}},{I^{e - 1}},{G^{e - 1}},{J^{e - 1}}] - t_{e - 2}^{1 - \varsigma }{\wp _2}[{t_{e - 2}},{\beta ^{e - 2}},{I^{e - 2}},{G^{e - 2}},{J^{e - 2}}]\} \times {\Xi _1}\\ &+ \frac{{\xi {{(\Delta t)}^\xi }}}{{2AB(\xi )\Gamma (\xi + 3)}}\sum \limits _{e = 2}^b \{ t_e^{1 - \varsigma }{\wp _2}[{t_e},{\beta ^e},{I^e},{G^e},{J^e}] - 2t_{e - 1}^{1 - \varsigma }{\wp _2}[{t_{e - 1}},{\beta ^{e - 1}},{I^{e - 1}},{G^{e - 1}},{J^{e - 1}}]\\ &+ t_{e - 2}^{1 - \varsigma }{\wp _2}[{t_{e - 2}},{\beta ^{e - 2}},{I^{e - 2}},{G^{e - 2}},{J^{e - 2}}]\} \times {\Xi _2}, \end{aligned} \end{aligned}$$
(53)
$$\begin{aligned} \begin{aligned} G(b + 1)&= G(0) + \frac{{1 - \xi }}{{AB(\xi )}}t_b^{1 - \varsigma }{\wp _3}[{t_b},\beta ({t_b}),I({t_b}),G({t_b}),J({t_b})]\\ &+ \frac{{\xi {{(\Delta t)}^\xi }}}{{AB(\xi )\Gamma (\xi + 1)}}\sum \limits _{e = 2}^b {{\wp _3}} [{t_{e - 2}},{\beta ^{e - 2}},{I^{e - 2}},{G^{e - 2}},{J^{e - 2}}]t_{e - 2}^{1 - \varsigma } \times \Xi \\ &+ \frac{{\xi {{(\Delta t)}^\xi }}}{{AB(\xi )\Gamma (\xi + 2)}}\sum \limits _{e = 2}^b \{ t_{e - 1}^{1 - \varsigma }{\wp _3}[{t_{e - 1}},{\beta ^{e - 1}},{I^{e - 1}},{G^{e - 1}},{J^{e - 1}}] - t_{e - 2}^{1 - \varsigma }{\wp _3}[{t_{e - 2}},{\beta ^{e - 2}},{I^{e - 2}},{G^{e - 2}},{J^{e - 2}}]\} \times {\Xi _1}\\ &+ \frac{{\xi {{(\Delta t)}^\xi }}}{{2AB(\xi )\Gamma (\xi + 3)}}\sum \limits _{i = 2}^b \{ t_e^{1 - \varsigma }{\wp _3}[{t_e},{\beta ^e},{I^e},{G^e},{J^e}] - 2t_{e - 1}^{1 - \varsigma }{\wp _3}[{t_{e - 1}},{\beta ^{e - 1}},{I^{e - 1}},{G^{e - 1}},{J^{e - 1}}]\\ &+ t_{e - 2}^{1 - \varsigma }{\wp _3}[{t_{e - 2}},{\beta ^{e - 2}},{I^{e - 2}},{G^{e - 2}},{J^{e - 2}}]\} \times {\Xi _2}, \end{aligned} \end{aligned}$$
(54)
$$\begin{aligned} \begin{aligned} J(b + 1)&= J(0) + \frac{{1 - \xi }}{{AB(\xi )}}t_b^{1 - \varsigma }{\wp _4}[{t_b},\beta ({t_b}),I({t_b}),G({t_b}),J({t_b})]\\ &+ \frac{{\xi {{(\Delta t)}^\xi }}}{{AB(\xi )\Gamma (\xi + 1)}}\sum \limits _{e = 2}^b {{\wp _4}} [{t_{e - 2}},{\beta ^{e - 2}},{I^{e - 2}},{G^{e - 2}},{J^{e - 2}}]t_{e - 2}^{1 - \varsigma } \times \Xi \\ &+ \frac{{\xi {{(\Delta t)}^\xi }}}{{AB(\xi )\Gamma (\xi + 2)}}\sum \limits _{e = 2}^b \{ t_{e - 1}^{1 - \varsigma }{\wp _4}[{t_{e - 1}},{\beta ^{e - 1}},{I^{e - 1}},{G^{e - 1}},{J^{e - 1}}] - t_{e - 2}^{1 - \varsigma }{\wp _4}[{t_{e - 2}},{\beta ^{e - 2}},{I^{e - 2}},{G^{e - 2}},{J^{e - 2}}]\} \times {\Xi _1}\\ &+ \frac{{\xi {{(\Delta t)}^\xi }}}{{2AB(\xi )\Gamma (\xi + 3)}}\sum \limits _{e = 2}^b \{ t_e^{1 - \varsigma }{\wp _3}[{t_e},{\beta ^e},{I^e},{G^e},{J^e}] - 2t_{e - 1}^{1 - \varsigma }{\wp _4}[{t_{e - 1}},{\beta ^{e - 1}},{I^{e - 1}},{G^{e - 1}},{J^{e - 1}}]\\ &+ t_{e - 2}^{1 - \varsigma }{\wp _4}[{t_{e - 2}},{\beta ^{e - 2}},{I^{e - 2}},{G^{e - 2}},{J^{e - 2}}]\} \times {\Xi _2}. \end{aligned} \end{aligned}$$
(55)

where

$$\begin{aligned} \Xi= & (b-e+1)^{\xi } -(b-e)^{\xi }\nonumber \\ \Xi _1= & (b-e+1)^{\xi } (b-e+3+2\xi )-(b-e)^{\xi }(b-e+3+3\xi ) \nonumber \\ \Xi _2= & (b-e+1)^{\xi }\Big [2(b-e)^2+(3\xi +10)(b-e)+2\xi ^2+9\xi +12\Big ]\nonumber \\ & -(b-e)^{\xi }\Big [2(b-e)^2+(5\xi +10)(b-e)+6\xi ^2+18\xi +12\Big ] \end{aligned}$$
(56)

The flowchart shown in Fig. 3 represents a complete scenario for the numerical technique used in this paper.

Fig. 3
figure 3

The proposed work flowchart.

Full size image

Results of proposed scheme

This section discusses the numerical simulation of the suggested approach for the diabetes model using the Mittage–Leffler kernel. The system parameter values and initial conditions are explained in31 where \(\Lambda =0.05\), \(\sigma =0.6\) l=1, e=0.4, b=0.07, p=0.3, \(\rho =0.6\), \(\nu =0.01\) and d=0.6. Using novel techniques, the efficacy of the theoretical results provided is verified. Figures 4, 5, 6, 7, 8, 9, 10 and 11 display solutions for each compartment at a fixed fractal dimension of the system with varying fractional order values. Matlab coding is used to find the numerical simulation for the fractional-order diabetes model. For every component of the model, distinct fractional order differentiation is used, as indicated by the fractional order \(\xi =(0.1,0.95,0.90,0.85)\) and fractal dimension \(\varsigma\). Higher values indicate larger memory effects, which are typical of biological systems with complicated, delayed reactions. These values also affect how the dynamics change with time.

In Fig. 4a increase \(\xi\) stabilizes beta cells quickly, reflecting efficient insulin secretion, while decrease fractional orders causes slower response, indicating dysfunction. \(\beta\) dynamics take longer to stabilize with increased oscillations in Fig. 4b, modeling weaker glucose regulation seen in diabetes. Figure 5a when increase \(\xi\) ensures rapid insulin response and decay, mimicking healthy regulation. Insulin stabilization Fig. 5b takes longer with prolonged peaks when decrease fractional orders, modeling impaired pancreatic response. Increase in fractional order ensures quick glucose stabilization Fig. 6a, reflecting normal uptake. Glucose stabilization Fig. 6b is slower with significant oscillations at low \(\xi\), suggesting weak compensatory mechanisms. In Fig. 7a, when increase in fractional order leads to efficient glucagon regulation, supporting glucose balance. Delayed stabilization and oscillatory behavior across \(\xi\) indicate inefficient glucagon control in Fig. 7b, worsening glucose imbalance. In Fig. 8a lead to faster stabilization and stronger beta cell dynamics, modeling efficient beta cell function. Figure 8b delay stabilization, reflecting reduced beta cell activity. Increase response speed but with minor oscillations Fig. 8c, simulating robust beta cell recovery. Figure 9a ensure faster stabilization, representing efficient insulin release and regulation. Figure 9b delayed stabilization and reduced peaks, simulating weakened insulin secretion. Figure 9c stabilize quickly, with minimal oscillations, reflecting strong regulatory feedback. Figure 10a stabilize quickly with minimal oscillations, simulating effective glucose uptake. Figure 10b delayed stabilization with fluctuations, representing inefficient glucose control. Figure 10c stabilize rapidly but show minor oscillations, modeling robust glucose uptake mechanisms. Figure 11a stabilize efficiently, ensuring effective glucose mobilization. Figure 11b lead to prolonged stabilization, reflecting weak regulation. Figure 11c stabilize quickly with minimal fluctuations, representing strong compensatory mechanisms. In Figs. 8, 9, 10 and 11) Higher \(\xi\) values enhance glucagon stabilization, ensuring efficient glucose regulation during fasting. Efficient glucagon stabilization reflects proper metabolic balance. Delays and oscillations mimic glucagon dysregulation, contributing to hyperglycemia. Variations in fractional orders significantly influenced the dynamics of the model. Higher orders \((\xi = 1)\) resulted in faster stabilization and minimal oscillations, which reflect efficient glucose-insulin regulation. In contrast, lower orders \((\xi = 0.95, 0.90, 0.85)\) led to delayed stabilization and pronounced oscillations, mimicking pathological states such as insulin resistance and glucagon dysregulation. These findings highlight the fractional-order model’s capability to capture the complexities of real-world diseases, offering insights into both normal and abnormal glucose homeostasis. Clinically, such models could enhance personalized diabetes management by tailoring interventions based on specific fractional parameters. Notable trends include faster stabilization and smoother dynamics observed at higher fractional orders, which indicate healthy regulation. Conversely, lower orders exhibited oscillations that mirrored diabetes-related dysfunctions. Unexpectedly, these lower orders also amplified delays, providing insights into glucagon dysregulation and hyperglycemia. Overall, the fractional-order model demonstrates superior performance compared to traditional integer-order models by capturing memory effects, leading to faster stabilization and reduced oscillations for \(\beta\)-cells, insulin, glucose, and glucagon. Integer-order models fail to replicate the delayed responses and oscillatory behavior typical of pathological states, which limits their clinical applicability. The flexibility of fractional orders allows for more accurate simulations of diabetes dynamics, enabling better predictive insights and the potential for personalized treatment strategies. The impact of fractional orders on proposed model, highlighting memory effects in diabetes regulation. Lower fractional orders result in smoother, delayed responses, capturing long-term biological dependencies that integer-order models fail to represent. These results align with real-world observations where glucose and insulin exhibit history-dependent behaviors, crucial in diabetes progression and treatment strategies. The fractional approach provides a more accurate representation of metabolic processes, offering insights for personalized diabetes management. The fractal-fractional approach with Mittag–Leffler kernels provides a more flexible and accurate model of diabetes dynamics by capturing memory and delayed responses inherent in biological systems, which traditional integer-order models miss. This approach allows for better simulation of long-term disease progression and varied therapeutic effects. Consequently, it offers insights into both chronic and acute phases of diabetes, enhancing predictive accuracy for management and intervention strategies.

Fig. 4
figure 4

Simulation of \(\beta\) with different fractional order.

Full size image
Fig. 5
figure 5

Simulation of I(t) with different fractional order.

Full size image
Fig. 6
figure 6

Simulation of G(t) with different fractional order.

Full size image
Fig. 7
figure 7

Simulation of J(t) with different fractional order.

Full size image
Fig. 8
figure 8

Simulation of \(\beta (t)\) at fractal 0.1 and different initial and fractional values.

Full size image
Fig. 9
figure 9

Simulation of I(t) at fractal 0.1 and different initial and fractional values.

Full size image
Fig. 10
figure 10

Simulation of G(t) at fractal 0.1 and different initial and fractional values.

Full size image
Fig. 11
figure 11

Simulation of J(t) at fractal 0.1 and different initial and fractional values.

Full size image

Conclusion

This study presents a novel mathematical model that combines the Mittag–Leffler kernel and fractional fractals. The diabetes model study aims to provide a sophisticated framework for understanding disease dynamics through advanced mathematical modeling techniques, ultimately contributing to better management strategies for diabetes and its complications. Chaos control and PID controllability in diabetes modeling help stabilize erratic glucose levels, ensuring a steady response to treatment. The use of the fractal fractional operator in diabetes modeling enhances the accuracy and reliability of predictions by better capturing the complexities of biological interactions, ensuring solution uniqueness, and providing flexible modeling capabilities that adapt to individual patient needs. A diabetes model using the fractal-fractional Mittag–Leffler function could contribute to future healthcare by enhancing predictive accuracy in glucose-insulin interactions, enabling more personalized treatment approaches. This model could also inform early intervention strategies by identifying irregular patterns that signal diabetes progression, ultimately supporting better long-term management and prevention. The fractional-order models for glucose-insulin monitoring arise from the complexity of Mittag–Leffler functions and fractal-fractional derivatives, which demand high computational resources. Implementing numerical methods like Newton polynomials adds to this expense due to iterative calculations. Efficient algorithms and parallel computing can mitigate these challenges. The study demonstrates the effectiveness of fractional-order modeling and PID control in stabilizing chaotic glucose-insulin dynamics, improving accuracy and system controllability. While the model offers significant insights, its reliance on idealized parameters may limit real-world applicability. Future research should focus on refining the model by incorporating real-world data to validate and adjust the parameters, enhancing its applicability in clinical settings. Additionally, exploring the integration of machine learning techniques could provide adaptive modeling capabilities, allowing for real-time adjustments based on individual patient responses. Investigating the impact of various lifestyle factors, such as diet and exercise, on glucose dynamics within the model could further improve its predictive power.