Comparing the Numerical Solution of Fractional Glucose–Insulin Systems Using Generalized Euler Method in Sense of Caputo, Caputo–Fabrizio and Atangana–Baleanu

: The purpose of this paper is to present a fractional nonlinear mathematical model with beta-cell kinetics and glucose–insulin feedback in order to describe changes in plasma glucose levels and insulin levels over time that may be associated with changes in beta-cell kinetics. We discuss the solution to the problem with respect to its existence, uniqueness, non-negativity, and boundedness. Using three different fractional derivative operators, the proposed model is examined. To approximate fractional-order systems, we use an efficient numerical Euler method in Caputo, Caputo–Fabrizio, and Atangana–Baleanu sense. Several asymptomatic behaviors are observed in the proposed models based on these three operators. These behaviors do not appear in integer-order derivative models. These behaviors are essential for understanding fractional-order systems dynamics. Our results provide insight into fractional-order systems dynamics. These operators analyze local and global stability and Hyers–Ulam stability. Furthermore, the numerical solutions for the proposed model are simulated using the three


Introduction
It is well known that the secretion of insulin and the biological effectiveness of insulin during periods of glucose loading are strongly influenced by (i) the number and functional capacity of pancreatic beta cells, as well as (ii) peripheral resistance to insulin action during glucose loading.Diabetes mellitus is characterized by a slowing of insulin release, combined with a diminished secretory response to insulin, attributable to peripheral resistance to insulin.These factors can contribute to the development of glucose intolerance, the main symptom of the disease.Several mathematical models have been proposed over the years as a means of explaining the relationship between insulin and glucose concentrations in plasma in response to a load of glucose in the body.
There have been numerous mathematical models proposed over the years, but the most recent is proposed by Ackerman et al. [1], which describes the glucose tolerance test using two autonomous first-order differential equations .
In this case, x(t) and y(t) are the insulin and glucose levels, respectively, at a given point in time.In x(t) and y(t), the concentration differences between fasting levels are represented.As shown in the equation above, a 1 y represents the net increase in insulin due to glucose exposure, a 2 x represents the rate of insulin removal, a 3 y represents the rate of glucose removal independent of insulin, and a 4 x represents the rate of glucose reduction due to insulin.For both x and y, c 1 and c 2 are positive constants corresponding to zero order rates of increase.As early as the 1970s, Molnar et al. [2] developed a model that incorporated a modified form of Equation ( 2) and another mathematical equation that was derived directly from Equations ( 1) and (2).A model system comprised of three nonlinear differential equations was developed by Bajaj and Rao [3,4] that used the following formulation to account for the kinetic properties of beta cells within the glucose-insulin feedback system: dz dt = S 5 y(T − z) + S 6 z(T − z) − S 7 z. ( As can be seen in Table 1, there are a number of parameters associated with the model (3).Glucose level decreases as a result of insulin production S 5 The rate at which β-cells divide in response to blood glucose S 6 β-cell growth rate due to dividing and non-dividing cells S 7 The rate at which β-cells decrease as a result of its current level c 1 An increase in x at a constant rate c 2 An increase in y at a constant rate T The total number of dividing and non-dividing cells Overall, a negative feedback loop between β-cells and the body's basal insulin and glucose concentrations has been established [5,6].The model proposed by Bajaj and Roa in [3] addressed the role of β-cells in regulating plasma insulin levels by Lenbury et al. [7].The rates S 1 , S 2 , c 1 are determined by the density of the beta cells.This concept of fractional calculus plays a significant role in many branches of mathematics in addition to being important when attempting to model real-world problems [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26].For further information regarding the use of fractional-order differential equations in biological studies, please refer to [27][28][29][30][31][32][33][34][35][36][37][38][39][40].There are a variety of fractional derivatives with nonsingular kernels proposed for describing nonlocal systems in a more accurate and practical way.These derivatives have been compared to fractional derivatives.Based on the type of function specified by exponentials and Mittag-Leflers [36,37], there are two very plausible possibilities among the existing definitions.
This study presents a fractional nonlinear mathematical model incorporating betacell kinetics and a glucose-insulin feedback system that allows us to describe plasma glucose and insulin level changes over time.The model solution is discussed in terms of its existence, uniqueness, non-negativity, and boundedness.To approximate fractional-order systems, we employ an efficient numerical Euler method in Caputo, Caputo-Fabrizio, and Atangana-Baleanu.We investigate steady-state stability at the local and global levels.In comparison to other existing methods, the proposed approach is easier to implement and more efficient due to its properties.Based on the simulation, it is confirmed that the modified model with three different fractional operators exhibits distinct asymptomatic behaviors that are not found in the integer-order model.Fractional derivatives account for the system's history, incorporating past states into current dynamics.This allows the model to capture long-term dependencies and the cumulative effect of glucose and insulin levels over time, which is not possible with integer-order models.The fractional-order model allows greater flexibility in fitting experimental data.By adjusting the order of the fractional derivatives, the model can better match the observed behavior of glucose and insulin levels.This captures subtleties integer-order models might miss.Fractional operators are inherently nonlocal, meaning they consider the entire history of a function rather than just its local behavior.This nonlocality provides a more accurate description of the glucose-insulin feedback mechanism, which is influenced by past fluctuations in glucose and insulin concentrations.
This article addresses the development of a comprehensive mathematical model that accurately represents the dynamic interactions between insulin and glucose levels in the body.It incorporates the kinetics of beta cells.Traditional integer-order models cannot capture biological systems' complex, nonlocal, and memory-dependent nature.This study aims to overcome these limitations by employing fractional calculus to develop such a nonlinear mathematical model.The challenge lies in ensuring model accuracy, stability, and practicality while incorporating fractional derivatives.In addition, it lies in providing a numerical method that is both efficient and easy to implement.This approach is expected to offer more realistic insights into the glucose-insulin feedback system and better match quantitative observations to existing models.

Preliminaries and Definitions
In [35], the classic Caputo fractional derivative is presented.However, the classical version may not work properly in nonlocal dynamics due to its singular kernel.Many new derivatives have been introduced to overcome this drawback, such as Caputo-Fabrizio (CF) [36] and Mittag-Leffler (AB) [37].In the CF derivative, the exponential function is used as the kernel, while in the AB derivative, the ML function is used.Here are the definitions of the classic Caputo operator, CF, and AB.

Definition 1 ([35]). For
The Caputo derivative is a generalization of the Riemann-Liouville fractional derivative and is used to describe processes that are locally continuous and globally smooth.Additionally, the fractional integral is given as follows: Model (3) can be written in sense of Caputo as: Definition 2 ([36]).For h ∈ H 1 (a, b) and 0 < ρ 1 < 1, the Caputo-Fabrizio derivative is defined as follows: . The corresponding fractional integral is CF and has the following properties: Model (3) can be written in sense of Caputo as: Definition 3 ([37]).For h ∈ H 1 (a, b) and 0 < ρ 1 < 1, the AB derivative of order ρ 1 is and E ρ 1 is the ML function.Also, the corresponding AB fractional integral is Model (3) can be written in sense of Caputo as: According to the following equation, the differentiation operators discussed above have a useful relationship to their corresponding fractional integrals Each of the Caputo, CF, and AB fractional derivatives has its own advantages and disadvantages.The Caputo fractional derivative is widely used for its simplicity and ease of interpretation, while the CF fractional derivative offers better numerical stability and accuracy.The AB fractional derivative provides a more flexible and generalized framework for fractional calculus, but it can be computationally demanding and complex to implement.

Euler Method Involving Caputo, Caputo-Fabrizio, and Atangana-Baleanu Operators
There are several numerical methods presented in this section for solving the models ( 5), ( 7) and ( 9) using the fractional Euler method (10) that employs the Caputo, Caputo-Fabrizio, and Atangana-Baleanu operators.As a result of applying one of the fractional integral operators ( 4), (6), or (8) to Equation (4), we obtain: t indicates the integral operator corresponding to the fractional derivatives of Caputo, CF and AB.The numerical scheme is designed based on a uniform grid on [0, b] with time steps of h = b−0 N .We will use z k to approximate z(t) at t = t k , where t k = 0 + kh and k = 0, 1 . . ., N. By discretizing Equation ( 10) using the Euler method [34], we determine the following formulas for Caputo, CF and AB. where According to Equation ( 11), the discrete form of Equation ( 5) is In terms of discrete form, Equation ( 7) can be expressed as follows using Equation ( 12): Equation ( 9) can be expressed in discrete form according to Equation ( 13) This paper enhances the understanding of discrete approaches for solving fractionalorder nonlinear differential equations, facilitating more informed decisions in modeling and simulation.The comparative analysis underscores the importance of choosing the right numerical method for solving fractional-order nonlinear differential equations.While the Euler method offers simplicity and speed, the predictor-corrector method provides superior accuracy and stability, making it suitable for more complex and precise applications.
Researchers can leverage this analysis to optimize their numerical simulations, achieving reliable and efficient solutions tailored to their specific needs.
By conducting a detailed comparative analysis of the Euler method and the predictorcorrector method [41], the paper can better highlight the contributions and practical implications of the proposed approach.This analysis will provide a comprehensive understanding of the strengths and limitations of different numerical methods for solving fractional-order systems, guiding researchers in choosing the most suitable method for their specific needs.As can be seen in the following representation, the fractional-order system ( 14) is as follows:

Properties of the Solutions
where * a D t is one of the Caputo, CF and AB fractional derivatives, and |N|.
Then the norm of the matrix M = [m ij [t]] is defined by It is now important to establish the existence and the uniqueness of the solutions of (2) in the region Ω × (0, T ′ ] with respect to Theorem 1.In the region Ω, it is sufficient to maintain that fractional-order systems' solutions (14) will exist and be unique if and only if W(0) = W 0 exists in the region Ω.
Proof.The following is the solution to system (4) Consequently, we are left with the following inequality as a result The following inequality is obtained after some calculations where Theorem 2. Considering c 2 ≥ S 4 η as a means to return the solutions to the system (14), which start at R 3 + , are non-negative.
Proof.Negative values are not considered to be physically meaningful, so they are not considered when attempting to solve problems in biology.Positive bounded solutions are the most realistic and biologically relevant solutions, so they are often sought out when trying to solve biological problems.As a result, this leads to * a D It is concluded from Lemmas 5 and 6 under [42] that the Caputo fractional model ( 14) solutions can be classified as semi-positive under the assumptions.
Theorem 3. Solutions to model (2), which begin in R 3 + , are uniformly bounded by Proof.In order to follow the [34] approach, the following method is adopted.In order to understand this equation, take into consideration that Therefore, for any δ > 0, the following is true: One can choose δ < min{S 4 + S 2 , S 7 + S 6 T}.
It follows from [24], Lemma 9, that we are able to Here, E q is the Mittag-Leffler function as used by Mittag-Leffler.Based on [42], with reference to Lemma 5 and Corollary 6, one can obtain the following result: This result leads to the uniform boundedness in the region W for the solution of the model (14).This means that any solution u of the model starting in W will be bounded for all t ≥ 0.Moreover, the solutions will remain bound for all subsequent times.

Stability of the Proposed Model Locally and Globally
E * = (x * , y * , z * ) provides the steady-state value for the system, where Here z * is the steady-state value satisfied by the cubic equation. where Consequently, the Jacobian matrix for system (14) at equilibrium can be calculated as follows: After the linearized system has been linearized, the characteristic equation P 1 (λ) is then solved as follows: where Lemma 1.A local asymptotic stability exists to the equilibrium point E * = (x * , y * , z * ).
Proof.The Routh-Hurwitz array for the Jacobian matrix J(E * ) has been implemented in the following manner: It is believed that Routh-Hurwitz conditions [43] hold if you verify that (a 1 a 2 −a 0 a 3 ) a 2 has the same sign as a 1 .It follows that the three eigenvalues have negative real parts.In view of the fact that a 0 > 0, a 1 > 0, a 2 > 0, a 3 > 0, then (a 1 a 2 −a 0 a 3 ) a 2 > 0 and a 1 a 2 > a 0 a 3 hold.This leads to the equilibrium point being an asymptotic stability point, as per Routh-Hurwitz, and thus fulfilling the Routh-Hurwitz stability standards.
Based on the results cited in [13], we can conclude that Hopf bifurcation occurs when i.e., 4∆ − tr 2 (J)

Sensitivity to Initial Values
This section discusses the sensitivity of fractional glucose-insulin systems to the initial values they are given.Time plots have been taken for three sets of initial conditions: (0.5, 0, 0.5), (2, 0, 2), (1, 0, 1).Sensitivity to initial values implies that changes in the initial values cause fluctuations that are initially negligible, but quickly build up, as shown in Figure 1.The sensitivity of fractional glucose-insulin systems to initial values is a critical aspect that influences their dynamic behavior.By examining the time plots for different initial conditions, we observe that even small changes can lead to substantial differences in system trajectories.This sensitivity highlights the importance of accurate initial condition assessment and tailored interventions in managing glucose-insulin dynamics effectively.

Proof.
x − x 1 = 0 I Hence, the results follow.

Numerical Simulation
Table 1 provides support for the analytical results as well as assessing the effectiveness of the control strategy for fractional-order values ρ 1 using the Euler method in the sense of the Caputo, Caputo-Fabrizio, and Atangana-Baleanu operators.Consider the following: Example 1.The discretized fractional order we are considering is as follows: depending on the initial conditions, E 0 = (x(0), y(0), z(0)) = (2, 0, 2).
According to these parameters, E 0 is asymptomatically stable.As predicted, this result leads to the unique disease-free equilibrium E 0 .The comparison of x(t), y(t), and z(t) numerical solutions derived using the Euler method in the sense of Caputo, Caputo-Fabrizio, and Atangana-Baleanu operators for different values of ρ 1 = 0.65, 0.75, 0.85, 0.95, and 1.Using Caputo Euler, Caputo-Fabrizio Euler, and Atangana-Baleanu Euler methods, Figures 2-4 compare the numerical solutions to x(t), y(t), and z(t).As a result of these parameters, E 0 appears to be asymptomatic, that is, the solution converges to E 0 .The result is an increase in insulin uptake and an increase in insulin-excitable tissue glucose uptake, which results in a reduction in glucose levels.According to Figures 2-4, different values of ρ 1 result in different glucose activities of insulin-excitable tissue glucose uptake and insulin dynamics.Various cases also allow us to plot blood glucose levels, insulin-excitable tissue glucose concentrations, and insulin absorption from excitable tissues over time.

Discussion
The proposed method in this manuscript offers a balanced approach, combining accuracy, stability, robustness, and computational efficiency.Its ability to handle fractional-order systems effectively makes it a superior choice for solving complex nonlinear differential equations, particularly in biological modeling.By highlighting these advantages, the manuscript not only underscores the method's contributions but also provides practical guidance for researchers in selecting the most appropriate numerical techniques for their specific applications.The glucose-insulin system, despite being linear, contains significant and interesting dynamics that are crucial for understanding metabolic regulation and diabetes.These dynamics make it a valuable example for testing numerical solution algorithms and demonstrating the practical applicability of the proposed method.The insights gained from this model can serve as a basis for extending the method to more complex and nonlinear systems.This will enhance its relevance and utility in various fields of biomedical research.By focusing on numerical solution algorithms and using the glucoseinsulin system as an illustrative example, the manuscript can provide a comprehensive evaluation of the proposed method.Additionally, extending the discussion to nonlinear systems, such as the Lorenz and Chen systems, will demonstrate the versatility and robustness of the approach.This will make it more relevant and valuable for researchers in various fields.By expanding the method to solve nonlinear chaotic systems like the Lorenz and Chen systems, the authors can demonstrate the versatility and robustness of their proposed numerical solution algorithms.This additional example will strengthen the paper by showing that the method is not limited to linear systems.It can effectively handle nonlinear and chaotic dynamics.Fractional-order systems, characterized by differential equations involving derivatives of non-integer orders, pose unique challenges and opportunities for numerical analysis.In addition to the Euler and predictor-corrector methods [46], numerous other algorithms have been developed to solve fractional-order differential equations.Each method has its strengths and is suited to different types of problems.This section provides an overview of some notable solution algorithms and their applications, highlighting how they compare to the methods discussed in this paper.Given the diversity of solution algorithms, it is essential to choose the method that most effectively fits the specific requirements of problem at hand.The comparison results highlight the trade-offs between different numerical methods.The Euler method is suitable for quick approximations and scenarios where computational resources are limited.However, for applications requiring high accuracy and stability, such as long-term simulations of the glucose-insulin system, the predictor-corrector method is preferable despite its higher computational cost.For real-time monitoring and control systems, where quick approximations are needed, the Euler method may suffice.However, for detailed studies and accurate biological systems modeling, the predictor-corrector method is recommended.These findings guide researchers in choosing the appropriate method based on their specific needs and constraints.By providing a detailed discussion of the comparison results, the authors can underscore the strengths and limitations of each numerical method, helping readers understand the practical implications of their findings and guiding future research in the field.Local stability ensures the system can withstand small disturbances or noise.For instance, in the glucose-insulin system, local stability means that minor fluctuations in glucose levels or insulin production will not lead to significant deviations from the desired state.Predictability is crucial for maintaining system performance and designing effective control strategies.Global stability implies robustness to large disturbances and varying initial conditions.For the glucose-insulin system, global stability ensures that even significant deviations in glucose levels or insulin production will eventually stabilize to the desired state.This robustness is essential for ensuring the long-term health and functionality of the biological system, making it resilient to unexpected changes and varying initial conditions.By analyzing both local and global stability, we gain a comprehensive understanding of system behavior.Local stability provides insights into the system's response to minor disturbances, while global stability guarantees robustness against significant deviations.Together, these stability analyses demonstrate the effectiveness and reliability of the proposed fractional-order glucose-insulin model in maintaining stable plasma glucose and insulin levels under various conditions.Figures 2-4 show the trajectories of the state variables x(t), y(t), and z(t) starting from initial conditions near E * .As seen in the figure, all state variables converge to their respective equilibrium values x * , y * , and z * over time.Numerical simulations support Lemma 1.The convergence of the state variables to the equilibrium point E * , as demonstrated in Figures 2-4, confirms the local asymptotic stability of E * .Minor discrepancies between analytical and numerical results can be attributed to numerical approximation errors and computation precision.Overall, the consistency between numerical and analytical results reinforces Lemma 1.The numerical simulations provide strong evidence supporting the system's Hyers-Ulam stability.The small deviations between the perturbed and unperturbed solutions, as shown in Figures 2-4, validate the system's robustness against small perturbations.Minor discrepancies between analytical and numerical results can be attributed to numerical approximation errors.Overall, the consistency between numerical and analytical results reinforces the validity of the Hyers-Ulam stability for this system.By following this structured approach, the authors can effectively demonstrate how their numerical results support analytical findings regarding Hyers-Ulam stability, thereby providing a comprehensive validation of their theoretical work.

Conclusions
This section discusses several efficient numerical methods.As an example, there is the fractional Euler method, which involves the Caputo, Caputo-Fabrizio, and Atangana-Baleanu operators.It is important to keep in mind that dynamic analysis is concerned with positivity, boundedness, and equilibrium.Numerous well-known results have been demonstrated to be instrumental in the linearization of the model over the past few decades, namely the Jacobian criterion and the Routh-Hurwitz criterion.Based on cell culture studies and the modified glucose-insulin feedback relationship, we approximately solve the fractional-order mathematical model for diabetes.Using the fractional method ( 14), life, individuality, non-negativeness, and boundedness can be solved.Comparing the Caputo fractional order (14) results with the integer model, the actual data calculated by the fractional method ( 14) are presented.Based on these comparisons, the fractional order model represents the scheme more accurately than the integer order model.The glucose and insulin levels were calculated according to experimental observations under different experimental conditions, as well as the level of protein malnutrition in the control group.With basal levels, both intravenous glucose tolerance and insulin sensitivity responses were obtained.It can be predicted from the above model that the beta cell function (represented by the parameter S 6 ) declines in protein malnutrition.Consequently, the association between S 6 and beta cell function in our model can be justified, since a reduction in S 6 implies an impaired insulin response to glucose.By increasing the glucose level, the beta cells would be more stressed, reducing their efficiency.Even though b is used as a denominator in our glucose equation, it would not result in singularities since it represents a stress reduction factor for the number of beta cells, which never becomes vanishingly small even under extreme physiological conditions.Moreover, we conducted a linear stability analysis of Equation ( 14) and found that stability does not change as the parameters S 6 are decreased.It suggests that diabetes is a continuous transition from a normal to pathological state.Consequently, the fractional model could be applied to other diabetes classes.Insulin is a major pathophysiological factor in this process, as is peripheral insulin sensitivity.

Conflicts of Interest:
The author declares no conflicts of interest.

Figure 1 .
Figure 1.A comparison between the numerical solution of x(t) using the generalized Euler method in sense of Caputo, Caputo-Fabrizio and Atangana-Baleanu at different values of ρ 1 .

Figure 2 .Figure 3 .
Figure 2. A comparison between the numerical solution of x(t) using a generalized Euler method in sense of Caputo, Caputo-Fabrizio and Atangana-Baleanu at different values of ρ 1 .

Figure 4 .
Figure 4.A comparison between the numerical solution of z(t) using a generalized Euler method in sense of Caputo, Caputo-Fabrizio and Atangana-Baleanu at different values of ρ 1 .

Funding:
The author extend his appreciation to the Deanship of Scientific Research at Northern Border University, Arar, KSA for funding this research work through the project number "NBU-FFR-2024-871-01".Data Availability Statement:No new data were created or analyzed in this study.

Table 1 .
Model parameters (3): description of the model parameters.