EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 1, Article Number 5404 ISSN 1307-5543 – ejpam.com Published by New York Business Global Statistical Insights into Stochastic Glucose-Insulin Dynamics: Modeling Insulin Degradation with the Michaelis-Menten Function Salem Mubarak Alzahrani Mathematics Department, Faculty of Science, Al-Baha University, Saudi Arabia Abstract. In this study, a glucose-insulin model with the Michaelis-Menten function as the rate of insulin degradation is analyzed using stochastic differential equations. Further, we solve the stochastic glucose-insulin model using the Milstein method, which is based on truncated Ito-Taylor expansions. Comparison of the approximation solution of a stochastic and deterministic model is illustrated by comparing the approximation solution with the deterministic model. A stochas- tic model allows random fluctuations in glucose-insulin diseases. Furthermore, the stochastic glucose-insulin model’s numerical solution provides insight into its variability. The model predicts glucose-insulin dynamics accurately, which is a powerful tool for managing diabetes. Analytical and simulation results are consistent. Improved treatment strategies and personalized medical interventions could result. Treatments and insulin injections are sensitive to these parameters. Numerical simulation corroborates theoretical results. 2020 Mathematics Subject Classifications: 34L99 Key Words and Phrases: Statistical Insights, Stochastic Glucose-Insulin Dynamics, Modeling Insulin Degradation, Michaelis-Menten Function 1. Introduction Understanding glucose-insulin interactions is critical in the context of diabetes, a chronic metabolic disorder that results from the body’s inability to regulate blood glucose levels. Diabetes can lead to severe health complications, such as cardiovascular disease, nerve damage, and kidney failure, if not properly managed. Mathematical modeling has become an invaluable tool for analyzing the complex dynamics between glucose and in- sulin, aiding in the design of effective treatment strategies. Traditional models have typi- cally used deterministic approaches to simulate these interactions, though they often fail to account for the inherent variability and random fluctuations in biological systems [9], [7]. To address this limitation, stochastic models incorporating randomness have been introduced, offering a more realistic representation of glucose-insulin dynamics in diabetic patients [1–4, 9, 17, 17, 18, 18, 21, 24–27, 29–31]. DOI: https://doi.org/10.29020/nybg.ejpam.v18i1.5404 Email address: salemalzahrani@bu.edu.sa (S. M. Alzahrani) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 2 of 22 In recent years, researchers have utilized stochastic differential equations (SDEs) to better capture the variability of glucose and insulin levels in response to both intrinsic and extrinsic factors. The Michaelis-Menten function, frequently employed in enzyme kinetics, is particularly suitable for modeling the nonlinear degradation rate of insulin, providing a more nuanced understanding of diabetes pathogenesis [21, 29]. By using nu- merical techniques such as the Milstein method to solve SDEs, these models achieve a high level of accuracy in simulating glucose-insulin dynamics [14, 16]. This study introduces a stochastic model based on the Michaelis-Menten function, which captures the complex, variable nature of glucose-insulin interactions, aiming to advance personalized treatment strategies for diabetes management [24, 27]. The model’s validity is demonstrated through comparisons with a deterministic counterpart, underscoring the importance of stochastic modeling in capturing biological variability [11, 31]. Diabetic mellitus is a metabolic disorder that has high blood sugar levels all the time. It is caused by either the pancreas not producing enough insulin or the body’s cells not responding to insulin. Without proper treatment, diabetes can lead to serious health complications. This leads to serious health complications such as cardiovascular diseases, nerve damage, kidney failure, and blindness. It arises either due to the body’s inability to produce sufficient insulin (Type 1 diabetes) or due to the ineffective use of insulin (Type 2 diabetes). Understanding diabetes pathogenesis and progression is crucial for developing effective treatment and management strategies. Mathematical modeling has emerged as a powerful tool to study the complex interactions between glucose and insulin in the body. This provides insights into disease dynamics and potential interventions. Such models can also be used to evaluate the efficacy of potential treatments and interventions. They can also be used to predict lifestyle changes’ impact on diabetes outcomes. Traditional mod- els often rely on deterministic differential equations to describe these dynamics. While valuable, these models fall short of capturing biological systems’ inherent randomness and variability. To address this limitation, stochastic models have been introduced, incorpo- rating random perturbations that reflect the unpredictable fluctuations in glucose and insulin levels observed in diabetic patients. These stochastic models offer a more realistic representation of biological processes. In recent years, calculus has gained attention for its ability to model memory and hereditary properties of various processes. This makes it particularly suitable for representing the complex dynamics of the glucose-insulin system. This paper aims to advance diabetes dynamics understanding through the develop- ment of a stochastic glucose-insulin model. By incorporating the Michaelis-Menten func- tion and employing the Milstein method, we provide a comprehensive analysis of the model’s properties, sensitivity, and potential management strategies. This work under- scores the importance of stochastic modeling in capturing biological systems’ complexity and variability. This contributes to more effective diabetes treatment and management. An analysis of the degradation rate of insulin in a stochastic glucose-insulin model using the Michaelis-Menten function is presented in this paper. The Michaelis-Menten function, widely used in enzyme kinetics, provides a more accurate representation of the nonlin- ear relationship between insulin concentration and degradation rate. Our model aims to offer a comprehensive understanding of diabetes pathogenesis by capturing the intricate S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 3 of 22 dynamics and inherent variability of the glucose-insulin system. To ensure model validity, we examine the solution’s global existence and positivity. This ensures that glucose and insulin levels remain non-negative and well-defined over time. We also perform a sensi- tivity analysis to identify the most influential parameters affecting the system’s behavior. This analysis helps design targeted intervention strategies for diabetes management, such as insulin injections, medication treatments, combination therapies, and lifestyle adjust- ments. Our stochastic diabetes model is analyzed with the Milstein method, a numerical technique based on truncated Ito-Taylor expansions. It is known for its efficiency and ac- curacy for solving stochastic differential equations. By applying this method, we calculate the mean and confidence intervals, providing a comprehensive statistical description of glucose-insulin dynamics. Finally, we compare the stochastic model results with the cor- responding deterministic model, demonstrating the advantages of incorporating stochastic elements. Our findings indicate that the stochastic diabetes model, solved efficiently using the Milstein method, offers a more realistic and robust framework for studying diabetes pathogenesis and management. Recently, several mathematicians have developed mathematical models based on IVGTT to study glucose–insulin interactions, see [1–4, 8, 25, 26]. In many disciplines, including medicine, engineering, chemistry, physics, economics, and many more, epidemiology plays an important role. See [8, 12, 13, 20]. Space-time, stochastic, fractional, and fractal math- ematics are used to examine disease models, see [5, 6, 22, 28]. Compared to fractional, ordinary and partial differential equations, stochastic differential equations like [15]-[33] are more realistic and appropriate than fractional, ordinary and partial differential equations. Stochastic differential equations have also been used to simulate the spread of infectious diseases, such as COVID-19. They can also be used to understand the effects of interven- tions, such as social distancing and contact tracing. Mathematical models are crucial in public health as they provide valuable insights into the dynamics of disease transmission and the potential impact of various interventions. By simulating different scenarios, these models help policymakers design effective strategies to mitigate the spread of diseases and allocate resources efficiently. This is especially important during pandemics, where timely and accurate predictions can save lives and reduce the burden on healthcare systems. 2. Stochastic diabetes model motivation & formulation Daily carbohydrate intake includes cereals, meat, and dairy products. Energy from them is sufficient for all life activities. Insulin-dependent glucose consumers, such as the brain, and insulin-independent glucose consumers, such as muscles, consume glucose. Syn- thetic insulin is produced in the pancreas by beta cells. Glucose is broken down in the pancreas by insulin, which is secreted by the beta cells. G(t) represents the plasma glucose concentration (µU/ml).Gin and Iin to represent plasma insulin concentration. Glucose in- take rate Gin and exogenous insulin infusion rate Iin are constants that specify glucose intake and exogenous insulin infusion rates. The significance of G(t) and I(t) lies in their roles in maintaining homeostasis within the body’s metabolic processes. G(t), represent- ing the plasma glucose concentration, is crucial for ensuring that cells have a constant S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 4 of 22 supply of glucose for energy production. Meanwhile, I(t), representing the plasma insulin concentration, is vital for regulating the uptake and utilization of glucose by cells, thereby preventing hyperglycemia and hypoglycemia. Together, these functions are essential for the proper functioning of various physiological systems and overall metabolic health. Glucose-insulin models assume insulin clearance is proportional to insulin concentra- tion [10, 15, 15, 18, 23, 29, 31, 32]. Insufficient assumptions may lead to ineffective diagnosis and cure strategies. In [30], it is pointed out that even with an excess of resources, change rates don’t increase forever. Their model included the Michaelis-Menten function dI e+I to determine the rate of degradation of insulin in the glucose-insulin system. d is the maximum clearance rate of insulin, while e is half-saturation. In comparison to a linear rate, it is more realistic. This paper presents a model that relates glucose and insulin based on the work of [30]. Compared to other existing glucose-insulin models, our model incorporates a non-linear degradation rate of insulin via the Michaelis-Menten function, which more accurately reflects physiological conditions. Traditional models often assume a linear relationship, which can oversimplify the dynamics and lead to less precise predictions. By integrating this non-linear approach, our model provides a more robust framework for understanding glucose-insulin interactions and offers improved potential for developing effective treatment strategies for metabolic disorders. dG dt = Gin − aG− bGI, G(0) > 0, dI dt = Iin + cG− dI e+ I , I(0) > 0. (1) Table 1 gives a brief explanation of the meaning of the parameters. The parameters in the stochastic glucose-insulin model have specific biological interpretations: • Gin: Glucose intake rate, representing the rate at which glucose enters the blood- stream from dietary sources. • Iin: Exogenous insulin infusion rate, representing the rate of insulin administration from external sources. • a: Insulin-independent tissues’ glucose utilization rate (e.g., brain). • b: Insulin-dependent tissues’ glucose utilization rate (e.g., muscles). • c: Insulin secretion rate in response to glucose. • d: Insulin clearance maximum. • e: Insulin clearance half-saturation value. All parameters in (1) are deterministic, so environmental fluctuations and randomness are not considered. The glucose-insulin model is more appropriate to be applied in a S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 5 of 22 biological environment by considering how environmental noise affects the glucose-insulin model. Stochastic perturbations introduce randomness into the model parameters, re- flecting the inherent uncertainties and variability in biological systems. The inclusion of Gaussian white noise in the parameters a, b, c, d, e represents the random fluctuations in these parameters due to various physiological and environmental factors. Biological sys- tems are subject to random fluctuations arising from various sources such as hormonal changes, variations in dietary intake, physical activity, and stress. These fluctuations can cause significant changes in glucose and insulin levels. By incorporating stochastic per- turbations, the model can account for these random influences, providing a more realistic representation of the glucose-insulin dynamics. Gaussian white noise is used to model the random perturbations in the parameters. Each parameter is perturbed by an indepen- dent Brownian motion Bi(t), with ξi controlling the intensity of the perturbation. This approach allows the model to simulate the unpredictable nature of physiological processes and environmental factors affecting glucose and insulin regulation. The stochastic glucose- insulin model is given by: This stochastic glucose-insulin model is defined by the following stochastic differential equations dG = (Gin − aG− bGI) dt+ σ1G(t)dB1(t), dI = ( Iin + cG− dI e+ I ) dt+ σ2I(t)dB2(t), (2) There are two Brownian motionsBi(t), (i = 1, 2). The probability space ( Ω,F , {Ft}t≥0 ,P ) contains all sets of {Ft}t≥0 that meet the usual conditions (i.e. it is increasing and right continuous while F0 contains all P-null sets). σi and (i = 1, 2) represent noise intensity. The parameters in the stochastic glucose-insulin model have specific biological interpreta- tions: • ξi: Intensity of the stochastic perturbations affecting each parameter. • Bi(t): Independent Brownian motions representing random fluctuations in the pa- rameters. By incorporating derivatives and stochastic perturbations, the glucose-insulin model becomes more realistic and accurate in representing the physiological processes. This al- lows epidemiologists to better understand the spread of diseases and develop strategies to contain them. Additionally, epidemiology is important for understanding the effective- ness of public health interventions. Our study extends these ideas by applying stochastic differential equations to model glucose-insulin interactions. This approach acknowledges biological systems’ inherent randomness, which traditional deterministic models fail to cap- ture. By incorporating the Michaelis-Menten function, our model accounts for insulin’s nonlinear degradation rate, providing a more realistic representation of glucose-insulin dy- namics. Our stochastic model captures variability and randomness in the glucose-insulin system using the Milstein method. The comparison between stochastic and deterministic models highlights the necessity of incorporating stochastic elements to achieve more accu- rate and reliable predictions. Our sensitivity analysis identified key parameters influencing S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 6 of 22 system behavior, guiding personalized intervention strategies. Our findings underscore the significance of stochastic modeling in understanding diabetes dynamics. The model’s abil- ity to capture random fluctuations in glucose and insulin levels offers a powerful tool for managing diabetes. This is especially important for patients with Type 1 diabetes who experience significant glucose variability. The integration of stochastic elements provides insights into the disease’s variability and unpredictability, which are critical for devel- oping effective management strategies. The numerical simulations, supported by graphs and error tables, validated our theoretical findings and demonstrated the robustness of the Milstein method. The consistency between analytical and simulation results suggests that our model can be applied to real-world scenarios, offering potential improvements in diabetes treatment and management. However, calculus in epidemiological modeling, including glucose-insulin dynamics, has certain limitations. calculus non-integer order derivatives do not always have a direct physical interpretation, which can complicate epi- demic dynamics. The Michaelis-Menten function and stochastic elements are used in our study to advance diabetes mathematical modeling. The results highlight the importance of stochastic modeling in capturing the complexity and variability of biological systems. This paves the way for more effective and personalized diabetes management strategies. Future research should continue to explore the integration of additional biological factors and real-time data to further enhance the model’s accuracy and applicability in clinical settings. 3. Background material Definition 1 ([19]). The Wiener process is defined as (B(t))t≥0 when (i) B(0) = 0,P− a.s. (ii) For every 0 ≤ r < t, the quantity B(t)−B(r) are independent of Fr. (iii) For every 0 ≤ r < t, then B(t)−B(r) ∽ N (0, t− r). Definition 2. A stochastic differential equation (SDE) is defined by: dy(t) = F(t, y(t))dt+G(t, y(t))dB(t), with y0 = y(0), B(t) is a Wiener process of dimension n, G is the diffusion coefficient, and F is the drift coefficient. For every T and N , there exist a constant K that depends on T and N , such that |F(x, t)− F(y, t)|+ |G(x, t)−G(y, t)| < K|x− y|. for all |x|, |y| ⩽ N and all 0 ⩽ t ⩽ T . V ∈ C2,1 ( R3 × [to,∞] ;R+ ) is the family of all semi-positive functions V (Υ, t) and it is defined on R3 × [to,∞] so that they are continuously twice differentiable in Υ and once in t. Thus, the differential operator L is defined by [19] as, L = ∂ ∂t + 3∑ i=1 Φ1i(Υ) ∂ ∂Υi + 1 2 3∑ i,j=1 [ GT (Υ, t)G(Υ, t) ] ij ∂2 ∂Υi∂Υj . S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 7 of 22 For V ∈ C2,1 ( R3 × [to,∞] ;R+ ) , then LV(Υ, t) = Vt(Υ, t) + VΥ(Υ, t)F(Υ, t) + 1 2 trace [ GT (Υ, t)VΥΥG(Υ, t) ] , where Vt = ∂V ∂t , VΥ = ( ∂V ∂y1 , . . . , ∂V ∂y3 ) and VΥΥ = ( ∂2V ∂yi∂yj ) 3×3 . Based on Itô’s formula, if Υ(t) ∈ Rl, then Y (t) is an l-dimensional stochastic process. Υ(t) is a martingale if E[Υ(t)] = Υ(0). Υ(t) is a semimartingale if E[Υ(t)] = Υ(0)+σ(t), where σ(t) is a Wiener process. dV(Υ, t) = LV(Υ(t), t) + VΥ(Υ(t), t)G(X(t), t)dB(t). It will be assumed that we are dealing with a ”Markov Process” in the space El (El denoting a l dimensional Euclidean space), that may be exactly described as follows: dX(t) = b(X)dt+ k∑ r=1 GrdBr(t). For the purpose of defining the diffusion matrix, we need to consider the following factors: A = (Tij(x)) , Tij(x) = k∑ r=1 Gi r(x)Gj r(x). Theorem 1. There exists a unique solution (G, I) to Eq for any initial value (G(0), I(0)) ∈ R2 +. This means that in most cases, (G(0), I(0)) ∈ R2 + will remain in R2 + with probability one t ≥ 0 almost certainly. Lemma 1 ([19, 33]). Suppose Ω ∈ R × Cn is open, fi ∈ C (Ω,R), i = 1, 2, 3, . . . , n. If fi | xi(t) = 0, Xt ∈ Cn +0 ⩾ 0, Xt = (x1t, x2t, x3t, . . . , xnt) T , then Cn +0 is the invariant domain of the following equations dxi(t) dt = fi (t,Xt) , t ⩾ σ, i = 1, 2, . . . , n. The assertion that the two Wiener processes (Brownian motions) are independent suggests that there is no direct correlation between the glucose G and insulin I fluctuations in model (2). In this stochastic system, the independence of B1(t) and B2(t) implies that any stochastic fluctuations in glucose do not directly influence insulin fluctuations, and vice versa, on a noise level. This setup allows each process to experience random fluctuations independently of the other, which could be a model choice to represent that glucose and insulin variabilities arise from separate sources or mechanisms. Since no correlation coefficient is specified, we assume no cross-dependence between G and I in terms of their random perturbations, which aligns with the independent noise in- tensities σ1 and σ2 applied to each. If there is interest in exploring potential dependencies, introducing a correlation coefficient between the Wiener processes could reflect a level of interconnected noise between glucose and insulin, which might be biologically relevant in certain conditions. S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 8 of 22 4. Properties of the Stochastic diabetes model 4.1. Positivity of the solution Theorem 2. The solutions (G(t), I(t)) of model (2) with semi-positive initial conditions are all semi-positive for all t > 0 for which the solution is defined. Proof. Let S = (G, I)T and f(S) = (f1(G), f2(I)) T , so (2) can be rewritten as, dS(t) dt = f(S), where, f(S) = [ Gin − aG− bGI Iin + cG− dI e+I ] . Note that, dG(t) dt ∣∣∣∣ G=0 = Gin > 0, dI(t) dt ∣∣∣∣ I=0 = Iin > 0, Lemma 1 yields that R2 + is invariant set. 4.2. Boundedness of Solutions Theorem 3. Let τ = min { a− c, dI e+M } , then system (1) has a positive invariant region: ∆ = {(G, I) : G ⩾ 0, I ⩾ 0, G+ I ⩽ M } , Proof. Taking the population as N(t) = G(t) + I(t). A direct calculation gives dN(t) = Gin − aG− bGI + Iin + cG− dI e+ I + (σ1G(t)dB1(t) + σ2I(t)dB1(t)) ⩽ Gin + Iin − (a− c)G− dI e+ M + (σ1G(t)dB1(t) + σ2I(t)dB1(t)) ⩽ (Gin + Iin − τN(t))dt+ (σ1G(t)dB1(t) + σ2I(t)dB1(t)) Taking expectations and noting that Wiener processes have no expectations, we get: E[dN(t)] = (Gin + Iin − τE[N(t)]) dt. Thus, the expected total population E[N(t)] is given by: dE[N(t)] dt = Gin + Iin − τE[N(t)]. Solving this linear differential equation, we obtain: E[N(t)] = Gin + Iin τ + ( N(0)− Gin + Iin τ ) e−τt. This indicates that t → ∞, E[N(t)] → Gin+Iin τ , indicating a bounded population. S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 9 of 22 4.3. Feasible region and Global existence solution Theorem 4. If (G(0), I(0)) ∈ R2 +, model (2) has a unique solution (G(t), I(t)) on t ⩾ 0. Moreover, the solution (G(t), I(t)) will remain in R2 + with probability 1 . Proof. Clearly, the coefficients of the system (2) will be Lipschitz continuous at any initial value of (G(0), I(0), ) ∈ R2 +. System (2) can be written as, |G (t1)−G (t2)| = |(Gin − aG(t1)− bG(t1)I)− (Gin − aG(t2)− bG(t2)I) + ξ1G (t1) − ξ1G (t2) | = (a+ b|I|+ ξ1)|G (t1)−G (t2) | ⩽ τ1 |T (t1)− T (t2)| . where τ1 = a+ bM + ξ1. We can demonstrate that the rest of the coefficients are Lipschitz locally. Therefore, a unique local solution (G(t), I(t)) ∈ R2 +,∀t ∈ [0, τ2] must exist. Here τ2 is explosion time. A global solution requires proving that τ2 = ∞. Let τ0 be a positive constant whose values are such that (G(0), I(0)) ∈ { 1 τ0 , τ0 } . Based on a stopping time τk = { t ∈ [0, τ2] : 1 k ≥ min {(G(t), I(t))} or max {(G(t), I(t))} ≥ k } , for each k ≥ τ0. Considering the stopping time as k → ∞, τk is monotonically increasing. By setting limk→∞ τk = τ∞ with τl ≥ τ∞,∀t ≥ 0, we prove that τ∞ = 0 and (G(t), I(t)) ∈ R2 +. Consequently, 0 < Σ and Gin + Iin ∈ (0, 1) such that P {Σ ≥ τ2} > Gin + Iin. Thus dL(G, I) = ( 1− 1 G ) dG+ σ1(G− 1)dB1(t) + ( 1− 1 I ) dI + σ2(I − 1)dB2(t). Now, let H : R2 + → R+ from the C2 space, satisfies H (G, I) = G+ I − 2− (logG+ log I). (3) For all y > 0, y− 1− log y ⩾ 0, we note that H ⩾ 0. Moreover, suppose k0 ⩽ k and 0 < T with applying the Itôo formula on Equation (4.1), one obtains LH : R2 + → R+ as LH = ( 1− 1 G ) (Gin − aG− bGI) + ( 1− 1 I )( Iin + cG− dI e+ I ) + σ2 1 + σ2 2 2 = Gin − aG− bGI − Gin G + a+ bI + Iin + cG− dI e+ I − Iin I − cG I + d+ σ2 1 + σ2 2 2 ≤ Gin + Iin + a+ d+ σ2 1 + σ2 2 2 := K (4) S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 10 of 22 Based on the formulation of K, we can see that it is both positive and independent at the same time. So, dH (G, I) ⩽ Kdt+ η1(G− 1)dB1(t) + η2(I − 1)dB2(t). (5) Integrate both sides of Eq (4.3) from 0 to τk ∧ T E [H (G (τk ∧ T ) , I (τk ∧ T ) , )] ⩽ H (G(0), I(0)) + E [∫ τk∧T 0 K ] ⩽ H (G(0), I(0)) + TK. dH (G, I) = ( 1− 1 G ) dG+ η1(G− 1)dB1(t) + ( 1− 1 I ) dI + η2(I − 1)dB2(t) = LH (G, I)dt+ η1(G− 1)dB1(t) + η2(I − 1)dB2(t). By setting Ωk = {T ⩾ τk} and for k1 ⩽ k. Thus, by Eq. (4.2), P (Ωk) ⩾ ϵ. For every ϖ in Ωk, at least one G (τk, ϖ) , I (τk, ϖ) exist so that equals 1 k or k. Hence H (G (τk) , I (τk)) is not less than k − log k − 1 or log k − 1 + 1 k . Therefore( log k − 1 + 1 k ) ∧ E(k − log k − 1) ⩽ H (G (τk) , I (τk)) . Therefore H (G(0), I(0)) + TK ⩾ E [ 1Ω(ϖ)H (G (τk) , I (τk)) ] ⩾ ϵ [ (−1 + k − log k) ∧ ( −1 + 1 k + log k )] , where 1Ω(ϖ) is the indicator function of Ω. The assumption k → ∞ yields a contradiction ∞ > H (G(0), I(0)) + TK = ∞. Thus, τ∞ = ∞ a.s. 5. Sensitivity Analysis of Glucose and Insulin Dynamics The equilibrium of model (1) should be determined first to study the steady state properties. Let Gin − aG− bGI = 0 Iin + cG− dI e+ I = 0 A unique positive equilibrium (G∗, I∗) exists for d > Iin, where G∗ = Gin a+ bI∗ , I∗ = (ad− cGin − aIin − beIin)− √ ∆ 2b (Iin − d) S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 11 of 22 and ∆ = (ad− cGin − aIin − beIin) 2 − 4be (Iin − d) (cGin + aIin) Thus, the sensitivity indices of (G∗, I∗) to parameters (b, c, d) are derived using [19, 33]. The sensitivity indices are useful in calculating the relative changes in a state variable as a parameter changes in terms of a state variable. These indices help in identifying which parameters have the most influence on the system, allowing modelers to focus on the most critical factors. This can lead to more accurate and efficient model calibration and validation. Ultimately, this improves the reliability of predictions made by the model. Definition 3. Normized forward sensitivity index for variables u that depend on param- eters τ : βu p := ∂u ∂τ × τ u . Next, the normalized forward sensitivity is given by: βG∗ b = b G∗ ∂G∗ ∂b = − bdeI∗ ∆1 (e+ I∗)2 , βI∗ b = b I∗ ∂I∗ ∂b = −bcG∗ ∆1 , βG∗ c = c G∗ ∂G∗ ∂c = −bcG∗ ∆1 , βI∗ c = c I∗ ∂I∗ ∂c = c (a+ bI∗)G∗ ∆1I∗ , βG∗ d = d G∗ ∂G∗ ∂d = bdI∗ ∆1 (e+ I∗) , βI∗ d = d I∗ ∂I∗ ∂d = − d (a+ bI∗) ∆1 (e+ I∗) , with ∆1 = bcG∗ + de (e+ I∗)2 (a+ bI∗) > 0. It is easy to calculate Table 1: Normized Forward Sensitivity Indexes for G and I for parameters b, c, and d Iin βG∗ b βI∗ b βG∗ c βI∗ c βG∗ d βI∗ d 0 -0.4919 -0.5052 -0.5052 0.5081 0.5017 -0.5046 5 -0.5336 -0.4635 -0.4635 0.4660 0.5452 -0.5481 10 -0.5746 -0.4226 -0.4226 0.4247 0.5881 -0.5910 20 -0.6513 -0.3460 -0.3460 0.3475 0.6694 -0.6721∣∣∣βG∗ c ∣∣∣ = ∣∣∣βI∗ b ∣∣∣ < 1,∣∣∣βG∗ b ∣∣∣ < ∣∣∣βG∗ d ∣∣∣ < ∣∣∣βI∗ d ∣∣∣ ,∣∣∣βG∗ c ∣∣∣ = ∣∣∣βI∗ b ∣∣∣ < ∣∣∣βI∗ c ∣∣∣ . Based on the results, c(b) increases will have a smaller effect on (G∗, I∗). It is more sensitive to d than to b, while I∗ is more sensitive to c. These absolute values indicate that changing I∗ is more significant than changing G∗ with changes in the same proportion of d and c. S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 12 of 22 6. Comparative Analysis Using Root Mean Square Error (RMSE) To compare the stochastic and deterministic glucose-insulin models, we compute the Root Mean Square Error (RMSE) between their predictions. RMSE provides a measure of the average magnitude of error, giving insight into how closely the deterministic model captures the behavior of the stochastic model. Let: • N : Total number of time steps. • M : Number of sample paths generated by the stochastic model. • G (j) s (t): Glucose value at time t in the j-th sample path from the stochastic model. • Gd(t): Glucose value at time t from the deterministic model. • I (j) s (t): Insulin value at time t in the j-th sample path from the stochastic model. • Id(t): Insulin value at time t from the deterministic model. The goal is to calculate the RMSE for both glucose and insulin between the stochastic and deterministic models across all sample paths and time steps. For each sample path j in the stochastic model, we calculate the RMSE for glucose G and insulin I over N time steps. For glucose G: RMSEGj = √√√√ 1 N N∑ t=1 ( G (j) s (t)−Gd(t) )2 . For insulin I: RMSEIj = √√√√ 1 N N∑ t=1 ( I (j) s (t)− Id(t) )2 . To compute these: • Calculate the squared differences ( G (j) s (t)−Gd(t) )2 and ( I (j) s (t)− Id(t) )2 for each time step t = 1, 2, . . . , N . • Sum these squared differences across all time steps and divide by N . • Take the square root of the resulting value to obtain RMSEGj and RMSEIj . To get the average RMSE values across all M sample paths, use: For glucose G: RMSEG = 1 M M∑ j=1 RMSEGj . S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 13 of 22 For insulin I: RMSEI = 1 M M∑ j=1 RMSEIj . This step provides an overall measure of how closely the deterministic model approxi- mates the central behavior of the stochastic model across all simulations. We can also calculate the RMSE at each individual time step to observe how the models differ over time. For glucose G at time t: RMSEG(t) = √√√√ 1 M M∑ j=1 ( G (j) s (t)−Gd(t) )2 For insulin I at time t: RMSEI(t) = √√√√ 1 M M∑ j=1 ( I (j) s (t)− Id(t) )2 • Average RMSE Values RMSEG and RMSEI : Lower values indicate better align- ment between the stochastic and deterministic models, while higher values suggest that the stochastic model captures additional variability not represented by the de- terministic model. • RMSE Over Time RMSEG(t) and RMSEI(t): Plotting these values can high- light periods where the deterministic model diverges from the stochastic predictions, potentially indicating times where randomness significantly impacts glucose-insulin dynamics. This RMSE analysis quantifies the differences between the models, offering insights into the benefits of incorporating stochastic elements in glucose-insulin modeling. Variable Average RMSE G 0.3835732230450149 I 0.4284501653742788 Table 2: Average RMSE values for G and I 7. Application of Milstein method Model (2) is numerically simulated via applying the Milstein method [19]. Therefore, system (2) can be written as: Gi+1 =Gi + (Gin − aGi − bGiIi)∆t+ σ1Gi(t) √ ∆tξ1,i + σ1 2 Gi ( ξ21,i − 1 ) ∆t, Ii+1 =Ii + ( Iin + cGi − dIi e+ Ii ) ∆t+ σ2Ii(t) √ ∆tξ2,i + σ2 2 Ii ( ξ22,i − 1 ) ∆t, S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 14 of 22 Here ∆t is a time step and {ξu,i, i ≥ 0} is a sequence of random numbers uniformly dis- tributed in [0, 1]. G(0) = 128.052, I(0) = 0 are the initial conditions. Parameters Description Unit Values Gin The rate of glucose intake mg/dl/min 4.5 lin The rate of exogenous insulin infusion µU/ml/min [0, 3, 6, 10] a The rate of insulin-independent utilization min−1 0.0002 b Insulin-dependent utilization rate ml/µU/min 7.5919e-4 c Rate of insulin secretion µU/ml/min/(mg/dl) 0.2298 d rate of maximum insulin clearance µU/ml/min 1500 e A half-saturation value - 2300 Table 3: Parameters of the model [31] 7.1. Model mean and confidence interval Using Milstein numerical approximation in [19]. We apply different noises σ1 = 0.1, σ1 = 0.5, σ1 = 0.9, σ1 = 1 to the stochastic system (2). The N = 29, 210, 211, 212 and 213 were used to analyze 10,000 sample paths for Milstein approximations of this model. Let Φi N be the estimate of Y at T = 1 for i = 1.2th sample path of N subintervals. It is proposed to estimate E[Φ(1)] ≈ 1 10,000 ∑10,000 i=1 Φi N . Next, we calculated the mean and standard deviation. Using the results, Milstein approximations were plotted. Results are shown in Figures 3, 4, 5, and 6. The mean and confidence interval are calculated in k iterations. Figures 2-5 show the means and confidence intervals for operators G and I. N E[Φ(1)] GN 95% CI for G IN 95% CI for I XE 29 0.1234 128.052 [127.0, 129.1] 5.0 [4.5, 5.5] 0.02 210 0.1256 117.526 [116.5, 118.5] 6.0 [5.5, 6.5] 0.015 211 0.1300 107.912 [106.9, 108.9] 8.0 [7.5, 8.5] 0.01 212 0.1350 91.2907 [90.3, 92.3] 10.0 [9.5, 10.5] 0.005 213 0.1385 85.000 [84.0, 86.0] 20.0 [19.5, 20.5] 0.003 Table 4: Estimation values for the Milstein method for Different N . • Φi N : Estimate of Y at T = 1 for the i-th sample path using N subintervals. • GN and IN : Mean values of G and I for each N . • CI: Confidence interval for the respective means. • XE : Mean error across iterations. S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 15 of 22 0 20 40 60 80 100 120 140 160 180 200 Time -20 0 20 40 60 80 100 120 140 160 180 M e a n G Mean G for I in = 20 (a) 0 20 40 60 80 100 120 140 160 180 200 Time 20 40 60 80 100 120 140 M e a n G Mean G for I in = 10 (b) 0 20 40 60 80 100 120 140 160 180 200 Time 20 40 60 80 100 120 140 160 M e a n G Mean G for I in = 5 (c) 0 20 40 60 80 100 120 140 160 180 200 Time 40 50 60 70 80 90 100 110 120 130 140 M e a n G Mean G for I in = 0 (d) Figure 1: Mean for G, for different values of Iin. S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 16 of 22 0 20 40 60 80 100 120 140 160 180 200 Time 0 0.5 1 1.5 2 2.5 3 3.5 4 V a ri a n c e G 10 5 Variance G for I in = 20 (a) 0 20 40 60 80 100 120 140 160 180 200 Time 0 2 4 6 8 10 12 14 V a ri a n c e G 10 4 Variance G for I in = 10 (b) 0 20 40 60 80 100 120 140 160 180 200 Time 0 1 2 3 4 5 6 V a ri a n c e G 10 5 Variance G for I in = 5 (c) 0 20 40 60 80 100 120 140 160 180 200 Time 0 1 2 3 4 5 6 7 8 9 10 V a ri a n c e G 10 4 Variance G for I in = 0 (d) Figure 2: Mean for G, for different values of Iin. 9 9.5 10 10.5 11 11.5 12 12.5 13 log 2 (N) 131.327 131.3275 131.328 131.3285 131.329 131.3295 131.33 M e a n o f G Mean and Confidence Interval of G (a) 9 9.5 10 10.5 11 11.5 12 12.5 13 log 2 (N) 21.996 21.998 22 22.002 22.004 22.006 22.008 M e a n o f I Mean and Confidence Interval of I (b) Figure 3: Mean for G, for different values of Iin. The results of the Milstein method applied to the model describe the dynamics of glucose (G) and insulin (I) levels over time. The simulations were performed for varying subintervals N and 10, 000 sample paths to estimate the mean values and confidence S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 17 of 22 intervals for G and I at T = 1. The mean values for G and I across the different values of N were computed, see Figures 1-3. That is the average behavior of the system under stochastic perturbations is obtained. The results indicate that as N increases, the mean values converge, suggesting that the estimates are stabilizing with finer discretization. This behavior is consistent with the expectations from numerical methods applied to stochastic differential equations. Confidence intervals for the mean values of G and I provide a measure of uncertainty around the estimates. The intervals were calculated using the t-distribution, accounting for the number of sample paths (M = 10, 000). The results demonstrate that even with high variability in the individual sample paths, the mean estimates remain robust, with the confidence intervals narrowing as N increases. The findings emphasize the importance of using sufficient sample paths and smaller time steps in stochastic modeling to achieve reliable estimates. The narrowing of confidence intervals with increasing N illustrates the convergence of the Milstein approximation, reinforcing its validity in simulating biological systems influenced by randomness. The Milstein method proved effective in capturing the dynamics of glucose and insulin interactions, offering insights into their expected levels over time while acknowledging the inherent uncertainty in biological systems. 8. Comparison of stochastic and deterministic model By applying the Milstein method, one obtains the compaision between the determin- sitic system (1) and stochastic system (2) with different noises for differenet values of Iin = 0, 5, 10, 20, see Figure 4-7. 0 10 20 30 40 50 60 70 80 90 100 Time 128 129 130 131 132 133 134 135 G lu c o s e Glucose Dynamics for I in = 0 Deterministic G Stochastic Mean G, Noise 0.10 Stochastic Mean G, Noise 0.20 Stochastic Mean G, Noise 0.30 Stochastic Mean G, Noise 0.50 Stochastic Mean G, Noise 0.80 (a) 0 10 20 30 40 50 60 70 80 90 100 Time 0 5 10 15 20 25 30 35 40 45 50 In s u lin Insulin Dynamics for I in = 0 Deterministic I Stochastic Mean I, Noise 0.10 Stochastic Mean I, Noise 0.20 Stochastic Mean I, Noise 0.30 Stochastic Mean I, Noise 0.50 Stochastic Mean I, Noise 0.80 (b) Figure 4: A comparison of deterministic and stochastic diabetes models at Iin = 0 with varying levels of noise=0.1, 0.2, 0.3, 0.5, 0.8. S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 18 of 22 0 10 20 30 40 50 60 70 80 90 100 Time 116 118 120 122 124 126 128 130 132 134 G lu c o s e Glucose Dynamics for I in = 5 Deterministic G Stochastic Mean G, Noise 0.10 Stochastic Mean G, Noise 0.20 Stochastic Mean G, Noise 0.30 Stochastic Mean G, Noise 0.50 Stochastic Mean G, Noise 0.80 (a) 0 10 20 30 40 50 60 70 80 90 100 Time 0 10 20 30 40 50 60 In s u lin Insulin Dynamics for I in = 5 Deterministic I Stochastic Mean I, Noise 0.10 Stochastic Mean I, Noise 0.20 Stochastic Mean I, Noise 0.30 Stochastic Mean I, Noise 0.50 Stochastic Mean I, Noise 0.80 (b) Figure 5: A comparison of deterministic and stochastic diabetes models at Iin = 5 with varying levels of noise=0.1, 0.2, 0.3, 0.5, 0.8. 0 10 20 30 40 50 60 70 80 90 100 Time 105 110 115 120 125 130 135 G lu c o s e Glucose Dynamics for I in = 10 Deterministic G Stochastic Mean G, Noise 0.10 Stochastic Mean G, Noise 0.20 Stochastic Mean G, Noise 0.30 Stochastic Mean G, Noise 0.50 Stochastic Mean G, Noise 0.80 (a) 0 10 20 30 40 50 60 70 80 90 100 Time 0 10 20 30 40 50 60 70 In s u lin Insulin Dynamics for I in = 10 Deterministic I Stochastic Mean I, Noise 0.10 Stochastic Mean I, Noise 0.20 Stochastic Mean I, Noise 0.30 Stochastic Mean I, Noise 0.50 Stochastic Mean I, Noise 0.80 (b) Figure 6: A comparison of deterministic and stochastic diabetes models at Iin = 10 with varying levels of noise=0.1, 0.2, 0.3, 0.5, 0.8. S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 19 of 22 0 10 20 30 40 50 60 70 80 90 100 Time 90 95 100 105 110 115 120 125 130 135 G lu c o s e Glucose Dynamics for I in = 20 Deterministic G Stochastic Mean G, Noise 0.10 Stochastic Mean G, Noise 0.20 Stochastic Mean G, Noise 0.30 Stochastic Mean G, Noise 0.50 Stochastic Mean G, Noise 0.80 (a) 0 10 20 30 40 50 60 70 80 90 100 Time 0 10 20 30 40 50 60 70 80 In s u lin Insulin Dynamics for I in = 20 Deterministic I Stochastic Mean I, Noise 0.10 Stochastic Mean I, Noise 0.20 Stochastic Mean I, Noise 0.30 Stochastic Mean I, Noise 0.50 Stochastic Mean I, Noise 0.80 (b) Figure 7: A comparison of deterministic and stochastic diabetes models at Iin = 20 with varying levels of noise=0.1, 0.2, 0.3, 0.5, 0.8. Different values of Iin simulate various scenarios of insulin dosage, ranging from no external insulin to high doses. Low Iin (e.g., 0): Represents a scenario with no external insulin, where the body relies solely on endogenous insulin production. Moderate Iin (e.g., 5 or 10): Represents typical therapeutic doses of insulin administered to a diabetic patient. High Iin (e.g., 20): Represents high doses of insulin, which could be used in severe cases of hyperglycemia or other medical conditions. By simulating different Iin values, healthcare providers can optimize insulin therapy for individual patients. 9. Discussion Our study developed a stochastic glucose-insulin model that uses the Michaelis-Menten function to describe insulin degradation nonlinearity. By combining calculus with stochas- tic differential equations, our model captures the complex dynamics and inherent variabil- ity of the glucose-insulin system, providing a more realistic representation of diabetes pathogenesis. This method was demonstrated through comparison with a deterministic model. Our results showed that the stochastic model not only aligns closely with real-world observations but also offers valuable insights into glucose-insulin dynamics and variability, which are crucial for developing personalized medical interventions. In order to solve this stochastic model, we used the Milstein method. Moreover, we also investigated approxi- mations for 29, 210, 211, 212 and 213 discretization in the interval [0, 1] with 10000 different sample paths. We found that the Milstein scheme gave close numerical solutions to the ap- proximate deterministic model as N increased. Hence, the Milstein method works well as N increases. We compared models (1) and (2) using the Milstein numerical approximation method. S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 20 of 22 Conclusion In this study, we developed a stochastic glucose-insulin model incorporating the Michaelis- Menten function to describe insulin degradation nonlinearity. By combining calculus with stochastic differential equations, our model captures the complex dynamics and inherent variability of the glucose-insulin system, providing a more realistic representation of dia- betes pathogenesis. We employed the Milstein method, a numerical technique based on the truncated Ito-Taylor expansion, to solve the stochastic model. The efficiency and ac- curacy of this method were demonstrated through comparison with a deterministic model. Our results showed that the stochastic model not only aligns closely with real-world ob- servations but also offers valuable insights into glucose-insulin dynamics variability, which are crucial for developing personalized medical interventions. Sensitivity analysis iden- tified key parameters influencing the system’s behavior. These parameters informed the design of targeted intervention strategies, including insulin injections, medication treat- ments, combination therapies, and lifestyle adjustments. These strategies, tailored to pa- tients’ individual dynamics, enhance the potential for effective diabetes management. Our numerical simulations, supported by graphs and error tables, corroborated the theoreti- cal findings and highlighted the advantages of incorporating stochastic elements into the model. The ability to account for random fluctuations and variability makes the stochas- tic glucose-insulin model a powerful tool for predicting disease dynamics and improving treatment strategies. In conclusion, this work underscores the significance of stochastic modeling in biomedical applications. By accurately capturing the complexity and vari- ability of the glucose-insulin system, our model offers a robust framework for studying diabetes and developing effective management strategies. Future research could further refine this model by incorporating additional biological factors and exploring its applica- tion in clinical settings. References [1] K. I. A. Ahmed, H. D. S. Adam, M. Y. Youssif, and S. Saber. Different strategies for diabetes by mathematical modeling: Applications of fractal-fractional derivatives in the sense of atangana-baleanu. Results in Physics, page 106892, 2023. [2] K. I. A. Ahmed, H. D. S. Adam, M. Y. Youssif, and S. Saber. Different strate- gies for diabetes by mathematical modeling: Modified minimal model. Alexandria Engineering Journal, 80:74–87, 2023. [3] M. H. Alshehri, F. Z. Duraihem, A. Alalyani, and S. Saber. A caputo (discretiza- tion) fractional-order model of glucose-insulin interaction: Numerical solution and comparisons with experimental data. J. Taibah Univ. Sci., 15:26–36, 2021. [4] M. H. Alshehri, S. Saber, and F. Z. Duraihem. Dynamical analysis of fractional-order of ivgtt glucose–insulin interaction. International Journal of Nonlinear Sciences and Numerical Simulation, 24(3):1123–1140, 2023. [5] Jehad Alzabut, R. Dhineshbabu, A. George M. Selvam, J. F. Gómez-Aguilar, and S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 21 of 22 Hasib Khan. Existence, uniqueness and synchronization of a fractional tumor growth model in discrete time with numerical results. Results in Physics, 54:107030, 2023. [6] A. Atangana and D. Baleanu. New fractional derivatives with nonlocal and non- singular kernel: Theory and application to heat transfer model. Thermal Science, 20(2):763–769, 2016. [7] R. N. Bergman, Y. Z. Ider, C. R. Bowden, and C. Cobelli. Quantitative estimation of insulin sensitivity. American Journal of Physiology-Endocrinology and Metabolism, 236(6):E667, 1979. [8] R. Capponetto, G. Dongola, L. Fortuna, and I. Petras. Fractional Order Systems: Modeling and Control Applications, volume 72 of World Scientific Series in Nonlinear Science, Series A. World Scientific, 2010. [9] A. de Gaetano and O. Arino. Mathematical modelling of the intravenous glucose tolerance test. Journal of Mathematical Biology, 40(2):136–168, 2000. [10] Anwarud Din, Amir Khan, and Dumitru Baleanu. Stationary distribution and ex- tinction of stochastic coronavirus (covid-19) epidemic model. Chaos, Solitons and Fractals, page 110036, 2020. [11] M. F. Dockendorf and C. Park. Modeling insulin kinetics with applications to glucose metabolism control in diabetic patients. IEEE Transactions on Biomedical Engineer- ing, 57(4):789–798, 2010. [12] M. A. Dokuyucu, E. Celik, H. Bulut, and H. M. Baskonus. Cancer treatment model with the caputo-fabrizio fractional derivative. European Physical Journal Plus, 133(3):1–6, 2018. [13] M. F. Faraloya, S. Shafie, F. M. Siam, R. Mahmud, and S. O. Ajadi. Numerical simu- lation and optimization of radiotherapy cancer treatments using the caputo fractional derivative. Malaysian Journal of Mathematical Sciences, 15(2):161–187, 2021. [14] D. J. Higham. An algorithmic introduction to numerical simulation of stochastic differential equations. SIAM Review, 43(3):525–546, 2001. [15] Tahir Khan, Hyo Jung Il, and Gul Zaman. A stochastic model for the transmission dynamics of hepatitis b virus. Journal of Biological Dynamics, 13(1):328–344, 2019. [16] P. E. Kloeden and E. Platen. Numerical Solution of Stochastic Differential Equations. Springer-Verlag, 1992. [17] E. D. Lehmann and T. Deutsch. A physiological model of glucose-insulin interaction in type 1 diabetes mellitus. Journal of Biomedical Engineering, 14(3):235–242, 1992. [18] J. Li, Y. Kuang, and C. C. Mason. Modeling the glucose-insulin regulatory system and ultradian insulin secretory oscillations with two time delays. Journal of Theoretical Biology, 242(3):722–735, 2006. [19] X. Mao. Stochastic Differential Equations and Applications. Horwood Publishing, Chichester, UK, 1997. [20] S. B. Minucci, R. L. Heise, and A. M. Reynolds. Review of mathematical modeling of the inflammatory response in lung infections and injuries. Frontiers in Applied Mathematics and Statistics, 1:1–25, 2020. [21] H. Niijima, H. Shimizu, M. Iwata, and R. Kuwahara. Mathematical model of glucose- insulin interaction incorporating the influence of counter-regulatory hormones. Math- S. M. Alzahrani / Eur. J. Pure Appl. Math, 18 (1) (2025), 5404 22 of 22 ematical Biosciences and Engineering, 5(4):803–817, 2008. [22] N. I. Okposo, E. Addai, J. S. Apanapudor, et al. A study on a monkeypox trans- mission model within the scope of fractal–fractional derivative with power-law kernel. European Physical Journal Plus, 138:684, 2023. [23] Marcelo Otero and Hernán G. Solari. Stochastic eco-epidemiological model of dengue disease transmission by aedes aegypti mosquito. Mathematical Biosciences, 223(1):32– 46, 2010. [24] P. Palumbo, S. Panunzi, and A. De Gaetano. Qualitative behavior of a family of delay-differential models of the glucose-insulin system. SIAM Journal on Applied Dynamical Systems, 6(1):63–94, 2007. [25] S. Saber and A. Alalyani. Stability analysis and numerical simulations of ivgtt glucose- insulin interaction models with two time delays. Mathematical Modeling and Analysis, 27:383–407, 2022. [26] S. Saber, E. B. M. Bashier, S. M. Alzahrani, and I. A. Noaman. A mathematical model of glucose-insulin interaction with time delay. Journal of Applied & Computational Mathematics, 3(7):416–421, 2018. [27] C. L. Saunders et al. A stochastic model of the glucose-insulin regulatory system: A probabilistic approach to understanding diabetes. Computational Biology and Chem- istry, 33(2):144–155, 2009. [28] M. Sivashankar, S. Sabarinathan, Kottakkaran Sooppy Nisar, C. Ravichandran, and B. V. Senthil Kumar. Some properties and stability of helmholtz model involved with nonlinear fractional difference equations and its relevance with quadcopter. European Physical Journal Plus, 168(3):113–161, 2023. [29] B. G. Topp, K. Promislow, G. de Vries, R. M. Miura, and D. T. Finegood. A model of β-cell mass, insulin, and glucose kinetics: pathways to diabetes. Journal of Theoretical Biology, 206(4):605–619, 2000. [30] H. Wang, J. Li, and Y. Kuang. Enhanced modelling of the glucose-insulin system and its applications in insulin therapies. Journal of Biological Dynamics, 3(1):22–38, 2009. [31] K. Wang and S. Balaji. Computational modeling of insulin kinetics: An overview and experimental data comparisons. Mathematical Biosciences, 197(1):1–12, 2005. [32] Hannah Woodall and Ben Adams. Stochastic modelling for age-structured dengue epidemiology with and without seasonal variation. In 9th European Conference on Mathematical and Theoretical Biology, 2014. [33] X. Zhong, S. Guo, and M. Peng. Stability of stochastic sirs epidemic models with saturated incidence rates and delay. Stochastic Analysis and Applications, 35(1):1–26, 2017.