Mathematics in Applied Sciences and Engineering https://ojs.lib.uwo.ca/mase Volume 4, Number 1, March 2023, pp.12-29 https://doi.org/10.5206/mase/15505 A NOVEL IVGTT MODEL INCLUDING INTERSTITIAL INSULIN JUN JIN, JIAXU LI∗, RUI XU, LEI YU, AND ZHEN JIN∗ Abstract. Minimal Model (MM) is the top-scoring model for assessing physiological characteristics to diagnose the potential or onset of type 2 diabetes mellitus (T2DM) through the intravenous glucose tolerance test (IVGTT) for the past four decades. Nevertheless it has been arguable that MM method either overestimates glucose effectiveness (GE) or underestimates insulin sensitivity (IS) in some cases by both biologists through in vivo experiments and mathematicians by analysis and/or simulations. We propose a novel model including the interstitial insulin according to physiology and adapted from the well accepted Sturis’ model for the glucose-insulin metabolic system suitable to the IVGTT setting. Our model consistently overcomes the aforementioned defects in a subgroup of subjects. In addition, the variable X for insulin action in MM might be appropriately interpreted as an increment of insulin in the interstitial space in response to the bolus stimulus, rather than being proportional to the interstitial insulin as believed. 1. Introduction Quantitatively assessing physiological characteristics, e.g., insulin sensitivity (IS) and glucose effec- tiveness (GE), is essential to diagnose the progression and/or onset of type 2 diabetes mellitus (T2DM) and develop drug for treatment. The gold standard for the assessment of IS is Hyperinsulinemic Eu- glycemic Glucose Clamp test originally proposed by DeFronzo et al.[11], which is direct and accurate, and widely accepted [32]. But it is very invasive, expensive, time consuming, and the subjects suffer. A much less invasive protocol as a result is the intravenous glucose tolerance test (IVGTT). The IVGTT data sampled at the time marks, for example, -10, -1, 0, 2, 3, 5, 7, 10, 13, 17, 21, 25, 30, 35, 40, 45, 50, 60, 70, 80, 90, 100, 120, 160 and 180 min after a rapid bolus intravenous glucose infusion, is used to estimate the parameter values of a differential equation model so that the aforementioned physiological characteristics can be assessed through these parameters. Many such models have been proposed at least as early as 1975 [28, 4, 5, 15, 33, 36, 27], but the top-scoring model widely used in experiments in both laboratory and clinic research is the minimal model (MM) developed by Bergman, Cobelli, and their colleagues a few years later [4, 5] for the past 40 years [2]. Commercial software include MIN- MOD, MINMOD Millennium [34, 7], and a software implemented in MLab [21] by Civilized Software Inc. (Silver Spring, MD 20906) are available. MM has been used by many biologists in their experiments [20]. The MM is given by Bergman et al. [4]:{ G′(t) = −[Sg + X(t)]G(t) + SgGb, X′(t) = −p2X(t) + p3[I(t) − Ib]+, (1.1) Received by the editors 15 November 2022; accepted 13 February 2023; published online 20 February 2023. 2020 Mathematics Subject Classification. Primary 54C40, 14E20; Secondary 46E25, 20C20. Key words and phrases. Diabetes, IVGTT, Minimal Model, remote insulin, time delay. ∗Corresponding authors. RX, LY and ZJ were supported in part by NSFC Grant 12271317, 61803242, and 12231012, respectively. ZJ was supported in part by Key R&D Project in Shanxi Province Grant 201803D31032. 12 https://ojs.lib.uwo.ca/mase https://dx.doi.org/https://doi.org/10.5206/mase/15505 A NOVEL IVGTT MODEL INCLUDING INTERSTITIAL INSULIN 13 with the initial condition G(0) = G0 and X(0) = 0, where G(t) and I(t) are respectively plasma glucose and insulin concentration, Sg is the GE index, the positive parameter p2 and p3 are used to calculate the IS index Si = p3/p2, and G0 is the initial value of the plasma glucose concentration value. [I(t)−Ib]+ = I(t)−Ib if I(t)−Ib ≥ 0; [I(t)−Ib]+ = 0 if I(t)−Ib < 0, where Gb and Ib denote the basal glucose and insulin level, respectively. The state variable X with initial condition X(0) = 0 was called insulin action in the original paper [4] and later was believed to be proportional to remote insulin (i.e. insulin in interstitial compartment) by, for example, Pacini and Bergman [34], Ader et al. [1], Bergman et al. [3] and, even recently, Bergman [2]. While along with the resounding success has been achieved by MM, several drawbacks and limitations in MM have been pointed out by both biologists and applied mathematicians [13, 39, 9, 15, 32, 37, 18, 14, 23]. Finegood and Tzur [13] pointed out that the estimation of Sg through MM may be not correct, later Caumo and Cobelli [8] elucidated that the inaccuracy is due to “the second hidden compartment”. Cobelli et al. [9] furthermore showed that MM overestimates Sg and underestimates Si. Patarrao et al. [37] believed that “many limitations of minimal model analysis stem from the fact that the model oversimplifies the physiology of glucose homeostasis”. Recently, Ha et al. [18], Fosam et al. [14] and Koh et al. [23] showcased that MM systematically underestimates Si in African American females comparing with non-Hispanic white females, which unveiled the paradox in predicting the risks of T2DM in different ethnic groups. In particular, Ha et al. [18] pointed out that the presence of a strong first insulin secretion phase leads to an underestimation of insulin sensitivity. After the creation of MM, deeper understanding for the action of interstitial insulin on blood sugar removal became much clearer. For example, an endothelial barrier delays the transcapillary transport of insulin from plasma to interstitial space [42] and capillary endothelium poses a barrier that delays the onset of muscle insulin action [49]. Factors that impede insulin access to muscle could contribute to increase insulin resistance [49]. In this paper, we propose a novel IVGTT model taking interstitial insulin as an explicit state variable based on the recent advanced physiology and the formulation from Sturis’ model in [46]. We demonstrate that the aforementioned drawbacks could be rectified in some subgroup of animal subjects. We organize this paper as follows. We present the model formulation in details in next section, followed by a section showing this model is well posed. Then in Section 4, we utilize Latin Hypercube Sampling method to generate pure random parameter values in the parameter space and then analyze the correlations of the profiles from MM, the new model and the nine available IVGTT data. We also elucidate that the new model would not produce ruinous dynamics when the model parameter values are within physiological ranges. We will discuss our findings in the last section. 2. Formulation of the novel IVGTT model adapted from a physiological model A well known model describes the general glucose-insulin metabolism was formulated by Sturis et al. [46], which contains three state variables for the plasma glucose (G), the plasma insulin (Ip) and the interstitial insulin (Ii) and three auxiliary variables mimicking time delays. Through this model, Sturis et al. [46] successfully elucidated and well accepted that the ultradian oscillation of insulin secretion in physiological settings is intrinsically caused by time delays and the transfer of insulin between the plasma compartment and interstitial compartment. A number of models in this area have been stemmed from Sturis’ model, particularly the metabolic feedback loop involving more organs [47], with explicit time delays [26, 24], and consequent attempts for algorithms of the artificial pancreas for type 1 diabetes (T1DM) by Wu et al. [50], Huang et al. [19], Kissler et al. [22] and Song et al. [44]. These, in addition to physiological observation by Prigeon et al. [39], allow us to confidently adapt the insulin in the interstitial compartment as an independent state variable as in Sturis’ model and the explicit time 14 J. JIN, J. LI, R. XU, L. YU, AND Z. JIN delay in the model in [26] to formulate a novel model suitable for the IVGTT environment for our aforementioned aims. Figure 1. Diagram of the model (2.1). Refer to the model diagram in Figure 1 for the glucose and insulin regulations, we denote G(t), Ip(t) and Ii(t) as the concentrations of the plasma glucose, plasma insulin and interstitial insulin at time t, respectively. According to the procedure of IVGTT, we assume that the glucose bolus is infused into the vein in the time interval [0, t0], where t0 = 2 or 3, and then our model is defined in [t0,∞). Insulin is secreted from pancreatic β-cells. Summarized by Straub and Sharp [45], glucose-stimulated biphasic insulin secretion in the IVGTT setting includes the KATP channel-dependent pathways, and KATP channel-independent pathways, respectively. The first phase of insulin secretion is resulted from the exocytosis of immediately releasable β-cell granules followed by the second phase of due to the KATP channel-independent pathways acting in synergy with the KATP channel-dependent pathway. In the same way as Bergman et al. [4], De Gaetano and Arino [15], Li et al. [25] and Panunzi et al. [36], we model the first phase quick insulin release by the initial condition Ip(t0) = Ip0 taken at the time mark t0. The subsequent long lasting and a larger amount insulin release in the second phase exhibits an explicit time delay in the same manner as in the physiological setting. Together we mimic the slower secretion through liver to plasma compartment by σf(G(t− τ)) with time delay τ > 0 and the maximum insulin secretion σ > 0. The insulin transfer relates to the biological action of insulin in the slowly equilibrating interstitial space, obeying a passive diffusion process determined by the difference between the concentration levels in the plasma and interstitial compartments, for which, we adapt the structure in Sturis’ model for insulin transfer between the plasma and interstitial space is described by the term E (Ip(t) − Ii(t)) with the diffusion transfer rate E > 0 between the two compartment. Different from the saturating responses of both insulin-dependent and insulin-independent glucose uptakes in the environment of the daily life [46, 26], under the three-hour short clinic setting of IVGTT, the response of β-cells to the bolus glucose injection into the vein is to release the insulin in the readily-releasable- pool abruptly followed by a large amount secretion [38, 10]. Thus, taking the same modeling approach A NOVEL IVGTT MODEL INCLUDING INTERSTITIAL INSULIN 15 as MM and most of the subsequent work [15, 16, 35, 27, 40, 41], we also assume the hepatic glucose production (HGP) to be a constant (b > 0), the insulin-independent glucose utilization to be a linearly dependent term SgG(t), and the insulin-dependent glucose removal in interstitial space to be a term of bi-linear dependence, SiG(t)Ii(t), where Sg stands for the GE index and Si stands for the IS index. The bolus glucose is infused in at t = 0 min and takes t0 minutes, typically t0 = 2 or 3, to complete. Thus our model is given as follows  G′(t) = b−SgG(t) −SiG(t)Ii(t), I′i(t) = E (Ip(t) − Ii(t)) −diiIi(t), I′p(t) = σf(G(t− τ)) −E (Ip(t) − Ii(t)) −dpiIp(t), (2.1) with initial condition G(t) = φ(t) > 0,φ ∈ C([t0 −τ,t0]), Ip(t0) = Ip0 > 0,Ii(t0) = Ii0 > 0, where Ip0 is the plasma insulin data at t0, and Ii0 is the corresponding interstitial insulin concentration. The initial function φ(t) = { Gb, for t ∈ [t0 − τ, 0), Gb + t −1 0 (G0 −Gb)t, for t ∈ [0, t0], (2.2) where G0 is the glucose data at the time mark t0 when the bolus glucose infusion is completed. Ac- cording to [17] and [26], the term f(G(t − τ)) = (G(t − τ))γ/(αγ + (G(t − τ))γ) models the insulin secretion stimulated by glucose G with time delay τ > 0, half saturation α > 0, and γ > 1. Positive constant parameters dpi and dii stand for the rates of insulin degradation in plasma and interstitial space, respectively. The model (2.1) is a generalization of the model in [27], which strengthens the stability properties for the models in [36] and [35]. It is easy to show that the model (2.1) assumes a unique positive equilibrium E0. In physiological observation, the dynamics of glucose, plasma insulin and interstitial insulin will return to their basal level after the test (180 min). We can therefore assume that the unique equilibrium point is at the basal level, that is, E0 = (Gb,Ipb,Iib) = (Gb,Ib,M1Ib), where M1 = E/(E + dii). This allows us to further express three model parameters in terms of other parameters, b = (Sg + SiM1Ib)Gb, σ = E(1 −M1 + dpi)Ib/f(Gb), dii = E(1 −M1)/M1, (2.3) which reduces the total number of parameters to be estimated from ten to seven. Notice that 0 < M1 = E/(E + dii) < 1. The relation Iib = M1Ib is in agreement with that the interstitial insulin concentration is a fraction of the concentration of the plasma insulin at equilibrium. The fraction is shown to be about M1 ≈ 0.60 according to [39, 42, 43]. So we assume throughout this paper that Iib < Ipb. (2.4) Remark 2.1. In MM, the initial condition of the insulin action X(0) = 0. This could hint X(t) might be interpreted as the increment of the interstitial insulin caused by the stimulus of abrupt glucose increase. We will discuss this in the last section in more details. In Model (2.1), the initial condition of the interstitial insulin concentration Ii0 is positive but a fraction of the plasma insulin at the basal state. Remark 2.2. Clearly the best way to verify the assumption of the state variable X in MM is through a carefully designed in-vivo or in human experiment by clinicians and biologists. However this had been 16 J. JIN, J. LI, R. XU, L. YU, AND Z. JIN overlooked in the past. We utilize the economic in-silico approach to perform the verification through a delay differential equation (DDE) model and its analysis. 3. Global stability of the equilibrium of the novel model We first ensure that the model (2.1) manifests ideal qualitative behaviours. By Lemma 1 in [27], the proof of the following proposition is straightforward. Proposition 3.1. All solutions of model (2.1) exist for all t > t0 and are strictly positive and bounded, that is, 0 < Ip(t) + Ii(t) ≤ MI . = max{Ip(t0) + Ii(t0), σ/ min{dpi,dii}} , for t > t0, and 0 < G(t) ≤ GM = max{G(t0), b/Sg} , for t > t0. Now we show that the equilibrium E0 is globally asymptotically stable by employing a suitable Lyapunov function. Theorem 3.2. If there exist positive constants A1,B1 and C1 such that the following assumptions (H1) and (H2) hold: (H1) τ < min { A1R + A1 K 2 − B1 L2 B1C , E + dpi − L2 C + D } , (H2) [ C1(E + dii) − B1Dτ + A1K 2 ][ E + dpi − τ(C + D) − L 2 ] − (B1 + C1) 2 4B1 E 2 > 0, where L = σ (γ + 1)2 4γα ( γ − 1 γ + 1 )γ−1 γ and R = Sg + SiIib,K = SiGM , C = RL/2,D = KL/2 and GM is the maximum value of G(t), then the equilibrium Eb of the system (2.1) is globally asymptotically stable. Proof. Assume that the condition (H1) and (H2) are satisfied. Let (G(t),Ip(t),Ii(t)) be any positive solution. Let V1(t) = A1 2 (G(t) −Gb)2 + B1 2 (Ip(t) − Ipb)2 + C1 2 (Ii(t) − Iib)2, and V2(t) = B1C ∫ t t−τ ∫ t z (G(s) −Gb)2dsdz + B1D ∫ t t−τ ∫ t z (Ii(s) − Iib)2dsdz. where A1,B1,C1 > 0. Define V (t) = V1(t) + V2(t). (3.1) We shall show that dV (t)/dt < 0 along any positive solution of the system (2.1) and thus the unique equilibrium point is globally asymptotically stable. We leave the details in Appendix A. � In Section 4, we will apply Theorem 3.2 to show the global stability in seven of the nine IVGTT experiments by determining the coefficients of the Liapunov function. Remark 3.1. If the parameters E = 0 and dii = 0, then the system (2.1) is reduced to the model (1) in [27], then the conditions (H1) and (H2) become A NOVEL IVGTT MODEL INCLUDING INTERSTITIAL INSULIN 17 (H1)′ τ < min { 2A1R + A1K − B1L B1RL , 2dpi − L L(R + K) } , (H2)′ (A1 − B1Lτ) [2dpi − τL(R + K) − L] > 0. Comparing with the conditions in Theorem 3 in [27], the two theorems might not imply each other. 4. Studies of the time courses X(t) and Ii(t) with IVGTT data 4.1. Data. We obtained nine individual data sets in IVGTT from literature [3, 34, 15, 36, 21] (see Table 1). We respectively fit MM (1.1) and the model (2.1) with these data to estimate the model parameters and then compare the profiles {X(tk)} from MM (1.1), the insulin profiles {Ip(tk)} and {Ii(tk)} from the model (2.1), and the IVGTT the data {I(tk)} for each subject and analyze their correlations, where {tk}, k = 1, 2, 3, · · · , are the time marks at which the data were sampled. We compare the indices of the IS and GE from MM (1.1) with that from the model (2.1) and then we observed noticeable improvement for a subgroup of subjects (dogs) in [3] labeled by Fig32A, Fig32B and Fig32C in this paper. 4.2. Method of generating the time courses. We first obtain the original parameter values and IS and GE indices Si and Sg of MM from the report for the subject 2 by MINMOD Millennium [6], and from [3] for the subject Fig32A, Fig32B and Fig32C. To obtain the estimated indices Si and Sg for the other subjects (6, 7, 8, 27, MLabEx) that are not available in the literature, we implemented MM (1.1) in Python, referenced from the details of the commercial IVGTT MM software implemented in MLab by Civilized Software Inc., Silver Spring, MD 20906 [21]. We implemented the model (2.1) in Python. Then, we employ Latin Hypercube Sampling (LHS) method [31] to estimate the parameter values for both models. All the parameter values are shown in Table 3 and the corresponding time courses are demonstrated in Figure 2 for the subject 2 and Figure 3 for the other subjects. LHS is an effective sampling technique for generating a nearly random sample of the parameter values from a multidimensional parameter space. The property of the stratified sampling ensures that LHS requires fewer samples to represent the real distribution than simple random sampling [29]. To ensure the randomness in parameter search by LHS, we search the parameter space for Sg,Si,p2 and p3 in the logarithmic manner for MM, while the searches of other parameter values are linear for our model (2.1). Among these samplings of parameter sets, we choose the parameter set that minimizes the deviation between the model profile and the IVGTT data characterized by the coefficient of determination (R- squared): R2 = 1 − ∑n i=1 (yexp(ti) −ysim(ti)) 2∑n i=1 (yexp(ti) −yexp) 2 where yexp(ti) and ysim(ti) are the IVGTT data and simulated profiles at time ti, respectively. To determine the parameter values of MM for the other subjects (6, 7, 8, 27 and MLabEx), the initial value of the glucose level G0 is also taken as a model parameter in addition to the three model parameters, Sg,p2 and p3. So the parameter space is given by ΘMM = {Sg,p2,p3,G0} , where the ranges of the parameter values are obtained from [7] listed in Table 2. The initial value of X is set to be zero. Linear interpolation of the sampled insulin data is performed and applied as the input I(t). The basal insulin level is taken as the insulin data at the time mark 180 min according to MLab’s implementation. Then we utilize LHS method to generate 1,000,000 independent sets of random parameter values. Then, the set of parameter values with the best fitting effect for the glucose data is taken as the parameter values from which the profiles G(t) and X(t) are determined. 18 J. JIN, J. LI, R. XU, L. YU, AND Z. JIN Table 1. Nine IVGTT data. The data of Subject 2 is from [34], Subjects 6, 7 and 8 are from [15], Subject 27 is from [36], Subject labeled by Fig32A, Fig32B and Fig32C are from [3], and Subject labeled by MLabEx is from [21]. S u b j 2 S u b j 6 S u b j 7 S u b j 8 S u b j 2 7 F ig 3 2 A F ig 3 2 B F ig 3 2 C M L a b E x m in G ( m g / d l) I ( p M ) G ( m g / d l) I ( p M ) G ( m g / d l) I ( p M ) G ( m g / d l) I ( p M ) G ( m g / d l) I ( p M ) G ( m g / d l) I ( p M ) G ( m g / d l) I ( p M ) G ( m g / d l) I ( p M ) G ( m g / d l) I ( p M ) 0 9 2 .0 0 1 1 .0 0 8 7 .7 4 6 7 .9 2 8 7 .2 1 3 8 .5 7 7 7 .9 9 5 7 .9 0 8 6 .4 7 4 4 .0 0 9 3 .3 0 2 2 .7 6 6 6 .5 9 5 .8 6 1 0 1 .4 3 2 8 .6 0 9 3 .0 0 3 .0 0 1 3 0 0 .0 1 5 0 .0 1 2 3 5 0 .0 0 2 6 .0 0 2 2 5 .4 7 4 1 3 .2 1 2 9 9 .3 7 1 7 9 .4 5 2 2 6 .4 2 1 0 3 1 .4 0 3 4 5 .9 0 1 0 3 6 .0 0 2 7 8 .5 2 9 9 .8 2 2 6 4 .9 2 7 7 .3 4 3 0 0 .0 0 1 1 2 .0 6 3 2 7 3 .5 1 1 1 0 .3 4 2 8 1 .3 8 1 2 5 .4 9 5 6 3 .0 0 1 5 5 .0 0 4 2 8 7 .0 0 1 3 0 .0 0 2 1 4 .1 5 4 1 0 .3 8 2 5 9 .9 6 1 0 3 .9 8 2 2 8 .9 3 9 1 5 .7 0 2 7 5 .6 4 1 0 6 7 .0 0 2 6 6 .3 5 8 4 .6 7 2 1 6 .2 3 7 8 .5 2 2 7 4 .2 3 8 4 .0 5 5 2 3 4 .8 5 8 5 .2 6 4 0 2 .0 0 1 7 3 .0 0 6 2 5 1 .0 0 8 5 .0 0 2 0 3 .7 7 3 0 5 .6 6 2 5 3 .2 5 9 9 .7 9 2 0 3 .7 7 7 5 9 .7 0 2 6 3 .0 3 9 1 4 .0 0 2 2 4 .1 1 7 6 .5 2 1 9 9 .0 5 5 5 .0 8 2 6 6 .3 6 9 6 .3 0 7 3 3 7 .0 0 1 2 2 .0 0 8 2 4 0 .0 0 5 1 .0 0 2 0 0 .0 0 2 8 6 .7 9 2 4 4 .0 3 9 3 .9 2 2 0 1 .2 6 7 7 2 .3 0 2 4 1 .4 1 4 1 5 .0 0 1 9 5 .4 7 7 3 .0 3 1 8 4 .0 1 5 6 .8 4 2 3 9 .9 1 9 1 .0 5 1 0 2 1 6 .0 0 4 9 .0 0 1 9 5 .2 8 2 3 4 .9 1 2 2 5 .5 8 1 0 4 .8 2 1 9 6 .2 3 6 4 6 .5 0 2 2 8 .8 0 4 5 5 .0 0 1 5 8 .2 3 7 9 .4 7 1 6 8 .2 6 6 5 .0 4 2 2 7 .7 4 9 5 .1 4 3 0 3 .0 0 8 3 .0 0 1 2 2 1 1 .0 0 4 5 .0 0 1 9 2 .4 5 3 1 7 .9 2 2 2 3 .9 0 7 7 .1 5 1 8 3 .6 5 6 6 9 .2 0 2 2 7 .9 0 4 0 4 .0 0 1 5 1 .0 7 7 1 .3 2 1 5 3 .9 4 5 0 .9 8 2 1 6 .2 9 1 1 2 .0 6 1 4 2 0 5 .0 0 4 1 .0 0 2 1 8 .8 9 2 1 6 .0 0 1 3 3 .1 7 4 6 .8 3 2 0 6 .2 7 9 2 .2 2 1 5 1 7 4 .5 3 2 7 8 .3 0 2 0 3 .7 7 8 8 .8 9 1 7 3 .5 8 5 1 3 .2 0 1 4 2 .4 8 6 0 .3 5 2 4 4 .0 0 6 3 .0 0 1 6 1 9 6 .0 0 3 5 .0 0 1 1 6 .7 1 3 1 .6 7 1 9 3 .3 8 1 0 3 .9 0 1 7 1 2 6 .7 3 4 3 .3 6 1 8 2 0 8 .9 8 3 4 4 .0 0 1 9 1 9 2 .0 0 3 0 .0 0 9 0 .9 3 2 2 .9 5 1 7 6 .9 2 9 9 .8 1 2 0 1 5 8 .4 9 2 3 8 .6 8 1 8 8 .6 8 9 5 .6 0 1 4 8 .4 3 5 0 8 .2 0 1 1 5 .2 7 3 3 .4 0 2 0 4 .0 0 4 4 .0 0 2 1 1 9 9 .9 7 2 8 2 .0 0 2 2 1 7 2 .0 0 3 0 .0 0 8 3 .7 7 2 0 .0 6 1 6 4 .7 4 8 9 .3 0 2 4 1 9 2 .7 7 2 3 2 .0 0 2 5 1 5 0 .0 0 2 5 0 .0 0 1 7 0 .2 3 7 9 .6 6 1 2 3 .2 7 4 4 0 .3 0 8 6 .6 3 1 8 .3 4 9 8 .0 9 2 5 .7 8 1 4 8 .9 9 8 5 .2 1 2 7 1 6 3 .0 0 2 7 .0 0 3 0 1 3 1 .1 3 2 3 3 .9 6 1 5 0 .9 4 9 7 .2 7 1 1 5 .7 2 3 2 7 .0 0 1 7 5 .6 5 2 9 4 .0 0 8 8 .7 8 2 4 .2 3 8 3 .0 5 2 0 .5 1 1 4 2 .5 0 6 7 .7 0 1 3 7 .0 0 2 8 .0 0 3 2 1 4 2 .0 0 3 0 .0 0 3 5 1 1 8 .8 7 2 0 3 .7 7 1 3 4 .1 7 8 6 .3 7 1 0 0 .6 3 2 8 6 .8 0 1 6 3 .9 4 1 9 3 .0 0 7 2 .3 2 1 4 .6 5 4 0 1 1 5 .0 9 1 5 3 .7 7 9 5 .6 0 2 2 6 .4 0 1 5 7 .6 4 2 2 7 .0 0 9 0 .9 3 1 9 .6 5 6 6 .5 9 1 3 .4 8 1 1 5 .2 4 3 9 .6 9 4 2 1 2 4 .0 0 2 2 .0 0 4 5 1 4 9 .5 3 2 1 0 .0 0 6 1 .5 8 8 .2 0 1 1 0 .0 0 2 0 .0 0 5 0 1 0 6 .6 0 1 6 9 .8 1 1 0 1 .4 7 4 4 .4 4 8 5 .5 3 1 6 6 .0 0 1 4 7 .7 3 1 8 8 .0 0 9 8 .8 1 1 4 .4 9 5 9 .4 3 8 .2 0 1 0 5 .8 4 2 6 .8 5 5 2 1 0 5 .0 0 1 5 .0 0 5 5 5 9 .4 3 8 .2 0 6 0 9 3 .4 0 1 1 5 .0 9 8 9 .7 3 2 4 .3 2 7 5 .4 7 1 4 8 .4 0 1 3 2 .4 1 1 1 6 .0 0 1 0 9 .5 5 1 2 .9 4 6 0 .8 6 8 .2 0 9 8 .9 0 8 .2 7 8 7 .0 0 1 0 .0 0 6 2 9 2 .0 0 1 5 .0 0 7 0 1 0 8 .9 9 1 9 4 .0 0 1 0 4 .5 4 1 1 .7 6 6 3 .7 2 6 .4 5 9 6 .7 0 1 3 .5 0 7 2 8 4 .0 0 1 1 .0 0 7 5 7 8 .0 0 4 .0 0 8 0 8 2 .0 8 1 1 1 .3 2 8 5 .5 3 3 3 .5 4 7 2 .9 6 1 1 8 .2 0 9 7 .2 8 1 5 4 .0 0 1 0 4 .5 4 1 2 .9 4 6 5 .1 6 4 .6 9 9 2 .3 4 2 8 .6 8 8 2 7 7 .0 0 1 0 .0 0 9 0 9 3 .6 8 9 5 .0 0 9 9 .5 2 1 9 .4 1 6 6 .5 9 5 .2 7 9 5 .8 6 6 .3 7 8 2 .0 0 2 .0 0 9 2 8 2 .0 0 8 .0 0 1 0 0 7 7 .3 6 5 3 .7 7 8 5 .5 3 2 9 .3 5 8 9 .1 8 7 2 .0 0 1 0 5 .9 7 1 1 .7 6 6 5 .8 7 4 .1 0 9 4 .3 7 5 .7 4 1 0 2 8 1 .0 0 1 1 .0 0 1 1 0 9 9 .5 2 1 5 .2 9 6 5 .1 6 3 .5 2 9 4 .3 2 1 5 .0 6 1 2 0 8 3 .0 2 4 6 .2 3 8 8 .0 5 3 7 .7 4 7 7 .9 9 6 7 .9 0 8 4 .6 7 5 0 .0 0 9 4 .5 1 1 1 .1 8 6 5 .8 7 4 .1 0 9 4 .9 7 2 1 .4 6 8 5 .0 0 6 .0 0 1 2 2 8 2 .0 0 7 .0 0 1 4 0 8 3 .0 2 5 8 .4 9 8 7 .2 1 3 1 .0 3 8 0 .5 0 4 2 .8 0 7 9 .2 7 3 8 .0 0 1 0 0 .2 4 1 8 .8 2 6 3 .7 2 6 .4 5 9 5 .5 7 6 .7 2 1 4 2 8 2 .0 0 8 .0 0 1 5 0 8 1 .0 0 1 .0 0 1 6 0 8 2 .0 8 6 4 .1 5 8 6 .3 7 3 3 .5 4 7 7 .9 9 6 0 .4 0 7 2 .0 6 3 6 .0 0 1 0 3 .1 0 1 1 .7 6 6 5 .1 6 5 .8 6 9 6 .1 8 2 0 .1 0 1 6 2 8 5 .0 0 8 .0 0 1 8 0 8 5 .8 5 5 5 .6 6 8 7 .2 1 4 6 .9 6 8 0 .5 0 5 7 .9 0 7 2 .0 6 3 3 .0 0 1 0 3 .0 0 1 1 .7 6 6 7 .3 0 5 .2 7 9 7 .4 9 8 .8 7 8 6 .0 0 1 .0 0 1 8 2 9 0 .0 0 7 .0 0 A NOVEL IVGTT MODEL INCLUDING INTERSTITIAL INSULIN 19 Table 2. Model parameter ranges. For MM (1.1), all the parameter ranges of MM are obtained from [7]. For Model (2.1), the ranges of Sg and Si are set slightly larger than that in [7] and the ranges of other parameter values are cited in the table. MM Value Range Model (2.1) Value Range Source Sg (1.2 × 10−3, 4.5 × 10−2) Sg (1 × 10−4, 9.9 × 10−2) [15, 7, 27] p2 (1.3 × 10−3, 2.0 × 10−1) Si (1 × 10−5, 9.9 × 10−3) [15, 7] p3 (5.4 × 10−7, 8.0 × 10−5) E (0.075, 0.25) [46] G0 (150, 400) τ (3, 25) [15, 27] dpi (0, 0.2) [26] α (100, 200) [48] γ (2, 5) [48] Table 3. Parameter values and IS and GE indices Si and Sg of MM and Model (2.1). Subj 2 Subj 6 Subj 7 Subj 8 Subj 27 Fig32A Fig32B Fig32C MLabEx MM by MLab: Si (×10−4) 7.46 0.63 4.02 0.382 0.27 9.10 5.40 1.00 4.21 Sg (×10−2) 0.95 1.40 1.85 2.36 1.33 7.70 4.70 5.20 1.86 G0 246.19 224.42 281.15 238.11 264.00 320.00 249.00 324.00 399.71 p2 (×10−2) 3.86 11.31 5.26 6.77 2.48 13.10 6.20 2.80 7.72 p3 (×10−5) 2.87 0.71 2.12 0.26 0.07 12.00 3.40 0.30 3.25 Model (2.1): Si (×10−4) 0.65 0.32 2.82 0.60 0.95 13.18 9.35 4.21 4.63 Sg (×10−2) 4.63 2.54 2.04 3.99 0.02 1.91 3.06 2.06 4.50 E 0.10 0.09 0.23 0.10 0.16 0.14 0.23 0.12 0.16 τ 6.06 6.34 9.30 8.72 15.70 3.33 5.31 6.51 21.49 dpi (×10−1) 1.58 1.54 1.33 1.39 0.89 1.64 1.48 0.99 0.40 α 102.20 135.95 106.00 118.34 124.00 177.00 164.27 181.93 140.93 γ 3.08 3.38 2.70 2.15 3.98 2.80 3.00 3.19 2.94 In the initial condition function (2.2) of the model (2.1), t0 = 2 or 3 (min) according to the experi- ments, and the basal levels of the glucose (Gb) and insulin (Ib) are determined by averaging the data at the time mark t = 0 and the last data point so that the measure errors for basal concentrations are reduced. Again we apply LHS method to produce 10,000 independent sets of random values in the parameter space Θ(2.1) = {Sg,Si,τ,E,dpi,α,γ} , and fit the glucose data for G and plasma insulin data for Ip to estimate the values of the seven parameters within the ranges given in Table 2. 4.3. Results. We take the subject 2 as an example to show the detailed comparison while summarize the comparisons of other subjects at the last. For the subject 2, applying Theorem 3.2 we find a set of coefficients of the Liapunov function V (t) =(G(t) −Gb)2 + (Ip(t) − Ipb)2 + (Ii(t) − Iib)2 20 J. JIN, J. LI, R. XU, L. YU, AND Z. JIN + 0.0017 ∫ t t−τ ∫ t z (G(s) −Gb)2dsdz + 6.8414 × 10−4 ∫ t t−τ ∫ t z (Ii(s) − Iib)2dsdz that ensures the global stability of the equilibrium E0, that is, both glucose and insulin concentrations will return their basal levels. We plot the IVGTT data, the profiles from MM, and the profiles from the model (2.1) in Figure 2. The dynamics of X(t) is multiplied by 4000 to fit into the scale of the window. Denote the values of the model profiles X, Ii, Ip using the parameter values in Table 3, and the IVGTT plasma insulin data at the sampling time marks tk by {X(tk)}, {Ii(tk)}, {Ip(tk)} and {I(tk)}, k = 1, 2, 3, · · · , respectively. Simple correlation analysis discloses that {Ip(tk)} and {I(tk)}, {Ii(tk)} and {I(tk)} are strongly correlated (0.989 and 0.909), but, on the contrast, {X(tk)} is not correlated with {Ii(tk)}, nor {I(tk)} with the correlation coefficients merely 0.520 and 0.219, respectively. 0 5 0 1 0 0 1 5 0 1 0 0 1 5 0 2 0 0 2 5 0 3 0 0 3 5 0 Glu cos e (m g/d l) t ( m i n ) G ( t ) G ( t ) g l u c o s e 0 3 0 6 0 9 0 1 2 0 1 5 0 I p ( t ) I i ( t ) X ( t ) p l a s m a i n s u l i n Ins uli n ( mU /m l) m o d e l ( 2 . 1 ) : d a t a : m o d e l M M : Figure 2. Profiles generated by MM and Model (2.1) for the subject 2. For six out of the other eight subjects, given the parameter values in Table 3, we are able to construct Liapunov functions in Theorem 3.2. The coefficients of Liapunov function defined by (3.1) are shown in Table 4 and hence the equilibrium E0 is globally asymptotically stable. We plot the IVGTT data and all model solution curves of these eight subjects in Figure 3 in smaller size of plots. Table 4. Coefficients of Liapunov functions in Theorem 3.2 for the subject 2, 7 and 8 in [27], the subjects in Fig. 32A, 32B, 32C in [3], and the example in [21]. Subj 7 Subj 8 Fig32A Fig32B Fig32C MLabEx A1/2 5 4 2 3 4 3 B1/2 3 1 2 1 3 2 C1/2 4 3 4 5 5 1 B1C 0.0132 0.0084 0.0056 0.0038 0.0078 7.6824 × 10−4 B1D 0.042 0.0028 0.068 0.028 0.0384 0.0044 The correlation analyses shown in Table 5 on the model profiles for the other eight subjects using the parameter values in Table 3 estimated in Section 4.2 consistently reveal that the plasma insulin profiles A NOVEL IVGTT MODEL INCLUDING INTERSTITIAL INSULIN 21 Table 5. Correlation coefficients between insulin concentrations in different compart- ments from MM (1.1) and the model (2.1). Subj 6 Subj 7 Subj 8 Subj 27 Fig32A Fig32B Fig32C MLabEx {Ip(tk), I(tk)} 0.958 0.939 0.929 0.961 0.967 0.960 0.975 0.984 {Ii(tk), I(tk)} 0.976 0.900 0.974 0.955 0.938 0.968 0.959 0.977 {Ip(tk), Ii(tk)} 0.986 0.984 0.977 0.984 0.991 0.991 0.990 0.996 {X(tk), I(tk)} 0.609 0.445 0.426 -0.112 0.576 0.526 0.146 0.492 {X(tk), Ii(tk)} 0.895 0.740 0.533 0.398 0.919 0.804 0.627 0.787 {Ip(tk)} from the model (2.1) are strongly correlated to the IVGTT plasma insulin data {I(tk)}. The interstitial insulin profile {Ii(tk)} by Model (2.1) is also strongly correlated to the plasma insulin data, although slightly weaker than the correlation between the plasma insulin {Ip(tk)} and the IVGTT data {I(tk)}. Without any surprise, the model profile {Ip(tk)} and {Ii(tk)} are strongly correlated as well. These analyses for the subject 2 and other eight subjects can be supported by the known physiological fact that interstitial insulin level is fractional of plasma insulin level, roughly about 60%, when at the basal state [39, 42, 43]. These reasonably indicate that the time course Ii(t) from the model (2.1) may indeed stand for the interstitial insulin concentration. On the other hand, neither of the nine profiles {X(tk)} of the auxiliary variable X in MM (1.1) show strong positive correlation to the sampled IVGTT plasma insulin data {I(tk)}. On the other hand, interestingly, the correlation coefficient between X and Ii are greater than 0.8 for Subject 6, Fig32A, Fig32B, three of the total nine subjects. But the three subjects inconsistently belong to different subgroups. The above analysis indicates that, even though X is not proportional to or consistently and strongly correlated to interstitial insulin, certain weak correlation between X and interstitial insulin still exists. In addition, observing that the initial condition X(0) = 0 and X(t) stabilizes at 0, perhaps, X(t) could be appropriately interpreted as a proportional increment of the interstitial insulin when the β-cells respond to the stimulus during IVGTT. Carefully examining the parameter values of the nine subjects in Table 3, it can be seen by comparison that the estimations by the model (2.1) is consistently overcome the aforementioned drawbacks of MM for the (dog) subject Fig32A, Fig32B and Fig32C [3]. For other (human) subjects, no consistent result is observed. Ha et al. [18] pointed out that the presence of a strong first insulin secretion phase could result in an underestimation of insulin sensitivity through a set of data generated by a putative model. Unfortunately we could not repeat their findings by our model (2.1). In a short summary, the above results manifest that a) the dynamics of Ip(t) and Ii(t) produced by the model (2.1) reflect the plasma and interstitial insulin levels during the IVGTT duration. b) in MM, the variable X(t) might be the increment of the interstitial insulin, which is not correlated to the interstitial insulin dynamics. Instead, it is more suitable to be considered as the proportional increment of the interstitial insulin from the basal level in responding to the bolus glucose stimulus. c) numerical studies for a subgroups (dog subjects) evidenced the potentiality of the model (2.1) to overcome the MM’s limitation – overestimating Sg and/or underestimating Si. 4.4. Physiological meaningful parameter values of E and τ would not destabilize the equi- librium. In this section, we consider the robustness of the model (2.1) in physiological applications. 22 J. JIN, J. LI, R. XU, L. YU, AND Z. JIN Figure 3. Profiles generated by MM and Model (2.1) for the other eight subjects. A NOVEL IVGTT MODEL INCLUDING INTERSTITIAL INSULIN 23 Under certain conditions, we have shown that the equilibrium E0 of our model (2.1) is globally asymp- totically stable. It is well known that larger time delays destabilize a system. We take the explicit time delay τ and the transfer rate E between the two insulin compartments as bifurcation parameters and investigate whether their changes would undermine the stability of E0. As results, omitting the routine but lengthy mathematical treatments, we obtained that, in addition to a natural physiological assumption Iib < Ipb (basal interstitial insulin level is lower than the basal plasma insulin level), a necessary condition for E0 to be unstable is Sg < (γ(1 −f(Gb)) − 1)SiIib. (4.1) For the parameter values of seven out of the nine subjects in Table 3 estimated from the IVGTT data, we successfully constructed Liapunov functions (refer to Table 4 for the coefficients) and thus the equilibrium point E0 is globally asymptotically stable. So we focus on the subject 27 and 6 now. The parameter values of both subjects do satisfy (4.1). However our intensive numerical simulations cannot detect a bifurcation value when E varies from 0 to large, but we did detect a Hopf bifurcation value for τ at τ0 ≈ 242.9 for subject 27 as shown in Figure 4, which is very much out of the physiological ranges and the instability would not happen. 230 235 240 245 250 255 60 70 80 90 100 20 30 40 50 60 70 tauG(t) I p (t ) Figure 4. Hopf bifurcation of Subject 27 when τ varies. Limit cycles bifurcated from a Hopf bifurcation when τ changes from 0 to large. The bifurcation point is at τ0 ≈ 242.9. 5. Discussions As the rapid increase of diabetic population in the world, research on the metabolic regulation in glucose and insulin becomes more and more pressing. Determining IS and GE are critical to find the progression pathways of T2DM and drug developments for T2DM. IVGTT is an appropriate protocol to determine these physiological characters less invasive and relatively accurate. However, the utility of the IVGTT for evaluating the essential physiological characteristics, insulin sensitivity and glucose effectiveness, has been challenged due to the minimal model [23]. A carefully formulated mathematical model for IVGTT setting is important to accomplish the assessment. One possible reason causing the 24 J. JIN, J. LI, R. XU, L. YU, AND Z. JIN drawbacks may be due to limited knowledge in biology regarding the interstitial insulin when the model was formulated, which resulted in a flaw in model structure. Till 2002, Straub and Sharp [45] pointed out that the glucose-stimulated biphasic insulin secretion in the IVGTT setting includes at lease two different pathways, the KATP channel-dependent pathways and KATP channel-independent pathways, which was not considered in the model formulations in MM. In this paper we formulate a reasonable model (2.1) to describe IVGTT based on physiological observations and the well accepted Sturis’ model for the glucose-insulin regulation. Our intensive numerical studies reveal that the parameter estimations of Si and Sg by our model show improvements of MM’s limitations for the dog subjects in the available data. Even though the number of the available data is limited, the consistency of the estimation in the entire subgroup may bring to light the right way of the model formulation. We analytically justified that the model (2.1) is well posed and the unique positive equilibrium point is globally asymptotically stable under certain conditions, which precludes the model (2.1) from the invalidity encountered by the extended MM by Pacini and Bergman [34] as pointed by De Gaetano and Arino [15]. Seven of the nine (≈ 78%) IVGTT data used in this paper do satisfy such conditions. Bifurcation analysis shows that the equilibrium point is stable when the values of the key parameter E and τ are within physiological ranges, which makes applications at ease. In addition, we also observed followings through our analyses. a) As known, in daily life (also called free living), the interstitial insulin and the plasma insulin keep at a balanced stable state [39, 42, 43]. In IVGTT, the balance is contravened for about 100 minutes due to the readily released insulin stored in β-cell granules into the peripheral triggered by the bolus glucose injection into the vein. b) The interstitial insulin level can be close to the level of the plasma insulin when the second phase insulin release occurs in the first 30 minute but hardly overpasses the plasma insulin level. The balance is gradually restored. c) When the plasma insulin exhibits the second peak caused by the second phase secretion, the interstitial insulin also acts like a bump but not a peak. In such case, the difference between the plasma insulin and interstitial insulin is larger. MM directly models the insulin action. When the insulin actions on adipocytes, skeletal muscle and liver in obesity are impaired, or decreased in the pathways of glucose uptake and metabolism, the glucose removal is eventually reduced [10]. Our analysis manifests that the variable X for insulin action in MM could be understood as a quantity proportional to the increment of the interstitial insulin, rather than proportional to nor strongly correlated to the interstitial insulin. It has been shown that hostile environmental factors promote the development of T2DM. Marunaka [30] found that the level of pH (potential of hydrogen) in the interstitial fluid is a critical factor re- sponsible for the occurrence of insulin resistance as pH level affects the binding affinity of insulin to its receptors. Even under mild metabolic disorder conditions with blood pH is kept constant within a normal range (7.35 − 7.45), the interstitial fluid pH would be lower than a normal level due to the less pH-buffering molecules in interstitial fluids than the powerful pH-buffering molecules in blood (e.g., hemoglobin and albumin), which leads to insulin resistance [30]. Quantifying the impact of pH on in- sulin resistance in detail by mathematical models could be significant in this area. Other environmental factors promote the development of T2DM can be found in the review by Dendup et al. [12]. The protocol IVGTT is less invasive, less expensive, but more accurate than other protocols such as oral glucose tolerant test (OGTT). If a reliable model can come up with for the estimation of the key physiological characteristics, it will be beneficial to both biologists and clinicians working in this area, and ultimately to the diabetic patients. A NOVEL IVGTT MODEL INCLUDING INTERSTITIAL INSULIN 25 Acknowledgements All the authors are grateful to the anonymous reviewers for their careful reviews and valuable com- ments. References 1. M. Ader, T. C. Ni, and R. N. Bergman, Glucose effectiveness assessed under dynamic and steady state conditions. comparability of uptake versus production components, J. Clin. Investig. 99 (1997), 1187–1199. 2. R. N. Bergman, Origins and history of the minimal model of glucose regulation, Front. Endocrinol. 11 (2021), 583016. 3. R. N. Bergman, D. T. Finegood, and M. Ader, Assessment of insulin sensitivity in vivo, Endocr. Rev. 6 (1985), 45–86. 4. R. N. Bergman, Y. Z. Ider, C. R. Bowden, and C. Cobelli, Quantitative estimation of insulin sensitivity, Am. J. Physiol. Endocrinol. Metab. 236 (1979), E667. 5. R. N. Bergman, L. S. Phillips, and C. Cobelli, Physiologic evaluation of factors controlling glucose tolerance in man: measurement of insulin sensitivity and beta-cell glucose sensitivity from the response to intravenous glucose, J. Clin. Investig. 68 (1981), 1456–1467. 6. Bergman-Laboratory, Minimal models for glucose and insulin kinetics, https://www.cedars-sinai.edu/research/ labs/bergman/resources.html;http://www.civilized.com/mlabexamples/glucose.htmld/, 2022, [accessed in Feb, 2022]. 7. R. C. Boston, D. Stefanovski, P. J. Moate, A. E. Sumner, R. M. Watanabe, and R. N. Bergman, Minmod millennium: a computer program to calculate glucose effectiveness and insulin sensitivity from the frequently sampled intravenous glucose tolerance test, Diabetes Technol. Ther. 5 (2003), 1003–1015. 8. A. Caumo and C. Cobelli, Minimal model estimate of glucose effectiveness: role of the minimal model volume and of the second hidden compartment, Am. J. Physiol. Endocrinol. Metab. 274 (1998), E573–E576. 9. C. Cobelli, A. Caumo, and M. Omenetto, Minimal model sgoverestimation and siunderestimation: improved accuracy by a bayesian two-compartment model, Am. J. Physiol. Endocrinol. Metab. 277 (1999), E481–E488. 10. M. P. Czech, Insulin action and resistance in obesity and type 2 diabetes, Nat. Med. 23 (2017), 804–814. 11. R. A. DeFronzo, J. D. Tobin, and R. Andres, Glucose clamp technique: a method for quantifying insulin secretion and resistance, Am. J. Physiol. Endocrinol. Metab. 237 (1979), E214. 12. T. Dendup, X.-Q. Feng, S. Clingan, and T. Astell-Burt, Environmental risk factors for developing type 2 diabetes mellitus: a systematic review, Int. J. Environ. Res. Public Health. 15 (2018), 78. 13. D.T. Finegood and D. Tzur, Reduced glucose effectiveness associated with reduced insulin release: an artifact of the minimal-model method, Am. J. Physiol. Endocrinol. Metab. 271 (1996), E485–E495. 14. A. Fosam, S. Yuditskaya, C. Sarcone, S. Grewal, H. Fan, and R. Muniyappa, Minimal model–derived insulin sensitivity index underestimates insulin sensitivity in black americans, Diabetes Care. 44 (2021), 2586–2588. 15. A. D. Gaetano and O. Arino, Mathematical modelling of the intravenous glucose tolerance test, J. Math. Biol. 40 (2000), 136–168. 16. A. D. Gaetano, D. D Martino, A. Germani, C. Manes, and P. Palumbo, Distributed-delay models of the glucose-insulin homeostasis and asymptotic state observation, IFAC Proceedings Volumes. 38 (2005), 1041–1046. 17. P. Gilon, M. A. Ravier, J. C. Jonas, and J. C. Henquin, Control mechanisms of the oscillations of insulin secretion in vitro and in vivo, Diabetes. 51 (2002), S144–S151. 18. J. Ha, R. Muniyappa, A. S. Sherman, and M. J. Quon, When minmod artifactually interprets strong insulin secretion as weak insulin action, Front. Physiol. 12 (2021), 508. 19. M. Huang, J. Li, X. Song, and H. Guo, Modeling impulsive injections of insulin: towards artificial pancreas, SIAM J. Appl. Math. 72 (2012), 1524–1548. 20. C. Izzi-Engbeaya, A. N. Comninos, S. A. Clarke, A. Jomard, L. Yang, et al., The effects of kisspeptin on β-cell function, serum metabolites and appetite in humans, Diabetes Obes. Metab. 20 (2018), 2800–2810. 21. Daniel R. Kerner, Minimal models for glucose and insulin kinetics, Civilized Software Inc., Silver Spring, MD 20906, https://www.civilized.com/pdffiles/glucose.pdf, 2021, [accessed in July, 2021]. 22. S. M. Kissler, C. Cichowitz, S. Sankaranarayanan, and D.M. Bortz, Determination of personalized diabetes treatment plans using a two-delay model, J. Theor. Biol. 359 (2014), 101–111. 23. H. C. E. Koh, B. W. Patterson, D. C. Reeds, and B. Mittendorfer, Insulin sensitivity and kinetics in african american and white people with obesity: Insights from different study protocols, Obesity. 30 (2022), 655—-665. 24. J. Li and Y. Kuang, Analysis of a model of the glucose-insulin regulatory system with two delays, SIAM J. Appl. Math. 67 (2007), 757–776. https://www.cedars-sinai.edu/research/labs/bergman/resources.html; http://www.civilized.com/mlabexamples/glucose.htmld/ https://www.cedars-sinai.edu/research/labs/bergman/resources.html; http://www.civilized.com/mlabexamples/glucose.htmld/ https://www.civilized.com/pdffiles/glucose.pdf 26 J. JIN, J. LI, R. XU, L. YU, AND Z. JIN 25. J. Li, Y. Kuang, and B. Li, Analysis of ivgtt glucose-insulin interaction models with time delay, Discrete Continuous Dyn. Syst. Ser. B. 1 (2001), 103. 26. J. Li, Y. Kuang, and C. C. Mason, Modeling the glucose–insulin regulatory system and ultradian insulin secretory oscillations with two explicit time delays, J. Theor. Biol. 242 (2006), 722–735. 27. J. Li, M. Wang, A. D. Gaetano, P. Palumbo, and S. Panunzi, The range of time delay and the global stability of the equilibrium for an ivgtt model, Math. Biosci. 235 (2012), 128–137. 28. V. Ličko and A. Silvers, Open-loop glucose-insulin control with threshold secretory mechanism: analysis of intravenous glucose tolerance tests in man, Math. Biosci. 27 (1975), 319–332. 29. J. Lin, R. Xu, and L. Li, Effect of leakage delay on hopf bifurcation in a fractional bam neural network, Int. J. Bifurcat. Chaos. 29 (2019), 1950077. 30. Y. Marunaka, Roles of interstitial fluid ph in diabetes mellitus: Glycolysis and mitochondrial function, World J. Diabetes. 6 (2015), 125. 31. M. D. McKay, R. J. Beckham, and W. J. Conover, A comparison of three methods for selecting values of input variables in the analysis of output from a computer code, Technometrics. 21 (1979), 21. 32. R. Muniyappa, S. Leen, H. Chen, and M. J. Quon, Current approaches for assessing insulin sensitivity and resistance in vivo: advantages, limitations, and appropriate usage, Am. J. Physiol. Endocrinol. Metab. 294 (2008), E15–E26. 33. R. V. Overgaard, J. E. Henriksen, and H. Madsen, Insights to the minimal model of insulin secretion through a mean-field beta cell model, J. Theor. Biol. 237 (2005), 382–389. 34. G. Pacini and R. N. Bergman, Minmod: a computer program to calculate insulin sensitivity and pancreatic responsivity from the frequently sampled intravenous glucose tolerance test, Comput. Methods Programs Biomed. 23 (1986), 113– 122. 35. P. Palumbo, S. Panunzi, and A. D. Gaetano, Qualitative behavior of a family of delay-differential models of the glucose-insulin system, Discrete Continuous Dyn. Syst. Ser. B. 7 (2007), 399–424. 36. S. Panunzi, Pasquale P. Palumbo, and A. D. Gaetano, A discrete single delay model for the intra-venous glucose tolerance test, Theor. Biol. Medical Model. 4 (2007), 1–16. 37. R. S. Patarrão, W. W. Lautt, and M. P. Macedo, Assessment of methods and indexes of insulin sensitivity, Revista Portuguesa de Endocrinologia, Diabetes e Metabolismo. 9 (2014), 65–73. 38. M. G. Pedersen and A. Sherman, Newcomer insulin secretory granules as a highly calcium-sensitive pool, PNAS. 106 (2009), 7432–7436. 39. R. L. Prigeon, M. E. Røder, D. Porte, and S. E. Kahn, The effect of insulin dose on the measurement of insulin sensitivity by the minimal model technique. evidence for saturable insulin transport in humans., J. Clin. Investig. 97 (1996), 501–507. 40. X. Shi, Y. Kuang, A. Makroglou, S. Mokshagundam, and J. Li, Oscillatory dynamics of an intravenous glucose tolerance test model with delay interval, Chaos. 27 (2017), 114324. 41. X. Shi, Q. Zheng, J. Yao, J. Li, and X. Zhou, Analysis of a stochastic ivgtt glucose-insulin model with time delay, Math. Biosci. Eng. 17 (2020), 2310–2329. 42. M. Sjostrand, S. Gudbjornsdottir, A. Holmang, L. Lonn, L. Strindberg, and P. Lonnroth, Delayed transcapillary transport of insulin to muscle interstitial fluid in obese subjects, Diabetes. 51 (2002), 2742–2748. 43. M. SjoStrand, S. GudbjoRnsdottir, L. Strindberg, and P. LoNnroth, Delayed transcapillary delivery of insulin to muscle interstitial fluid after oral glucose load in obese subjects, Diabetes. 54 (2005), 152–157. 44. X. Song, M. Huangn, and J. Li, Modeling impulsive insulin delivery in insulin pump with time delays, SIAM J. Appl. Math. 74 (2014), 1763–1785. 45. S. G. Straub and G. W. G. Sharp, Glucose-stimulated signaling pathways in biphasic insulin secretion, Diabetes Metab. Res. Rev. 18 (2002), 451–463. 46. J. Sturis, K. S. Polonsky, E. Mosekilde, and E. V. Cauter, Computer model for mechanisms underlying ultradian oscillations of insulin and glucose, Am. J. Physiol. Endocrinol. Metab. 260 (1991), E801–E809. 47. I. M. Tolić, E. Mosekilde, and J. Sturis, Modeling the insulin–glucose feedback system: the significance of pulsatile insulin secretion, J. Theor. Biol. 207 (2000), 361–375. 48. M. Wang, J. Li, G. E. Lim, and J. D. Johnson, Is dynamic autocrine insulin signaling possible? a mathematical model predicts picomolar concentrations of extracellular monomeric insulin within human pancreatic islets, PLoS One. 8 (2013), e64860. 49. I. M. Williams and D. H. Wasserman, Capillary endothelial insulin transport: The rate-limiting step for insulin- stimulated glucose uptake, Endocrinology. 163 (2022), bqab252. 50. Z. Wu, C. K. Chui, G. S. Hong, and S. Chang, Physiological analysis on oscillatory behavior of glucose–insulin regulation by model with delays, J. Theor. Biol. 280 (2011), 1–9. A NOVEL IVGTT MODEL INCLUDING INTERSTITIAL INSULIN 27 Appendix A. Details of the proof of Theorem 3.2 To calculate the derivation of V (t) = V1(t) + V2(t) along any positive solutions of system (2.1), we first have dV1(t) dt = A1(G(t) −Gb)[b−SgG(t) −SiG(t)Ii(t)] + B1(Ip(t) − Ipb) [σf(G(t− τ)) −E (Ip(t) − Ii(t)) −dpiIp(t)] + C1(Ii(t) − Iib) [E (Ip(t) − Ii(t)) −diiIi(t)] = A1(G(t) −Gb)[−Sg(G(t) −Gb) −SiG(t)(Ii(t) − Iib) −SiIib(G(t) −Gb)] + B1σ(Ip(t) − Ipb) (f(G(t− τ)) −f(Gb)) − B1(Ip(t) − Ipb) [(E + dpi) (Ip(t) − Ipb) −E(Ii(t) − Iib)] + C1(Ii(t) − Iib) [E(Ip(t) − Ipb) − (E + dii) (Ii(t) − Iib)] . (A.1) According to the mean value theorem, we have from (A.1) that dV1(t) dt = A1(G(t) −Gb)[−(Sg + SiIib)(G(t) −Gb) −SiG(t)(Ii(t) − Iib)] + B1σ(Ip(t) − Ipb)f′(G(ξ))(G(t− τ) −Gb) − B1(Ip(t) − Ipb) [(E + dpi) (Ip(t) − Ipb) −E(Ii(t) − Iib)] + C1(Ii(t) − Iib) [E(Ip(t) − Ipb) − (E + dii) (Ii(t) − Iib)] = − A1(Sg + SiIib)(G(t) −Gb)2 −A1SiG(t)(G(t) −Gb)(Ii(t) − Iib) − B1(E + dpi)(Ip(t) − Ipb)2 + (B1 + C1)E(Ip(t) − Ipb)(Ii(t) − Iib) − C1(E + dii)(Ii(t) − Iib)2 + B1σ(Ip(t) − Ipb)f′(G(ξ))[(G(t− τ) −G(t)) + (G(t) −Gb)], (A.2) where G(ξ) is between G(t−τ) and Gb. Because G is bounded and α and γ are positive constants, we get that f′(G(ξ)) is bounded. Similarly, setting R = Sg + SiIib,K = SiGM and L = σ (γ + 1)2 4γα ( γ − 1 γ + 1 ) γ−1 γ , where GM and L are the maximum values of G(t) and function of σf ′(G(ξ)) for G > 0, we obtain from (A.2) that dV1(t) dt ≤−A1R(G(t) −Gb)2 −A1K | (G(t) −Gb)(Ii(t) − Iib) | −B1(E + dpi)(Ip(t) − Ipb)2 + (B1 + C1)E(Ip(t) − Ipb)(Ii(t) − Iib) −C1(E + dii)(Ii(t) − Iib)2 + B1L | (Ip(t) − Ipb)(G(t− τ) −G(t)) | + B1L | (Ip(t) − Ipb)(G(t) −Gb) | . (A.3) Furthermore, since xy ≤ 1 2 (x2 + y2) for all x,y ≥ 0, we obtain that |(Ip(t) − Ipb)(G(t− τ) −G(t))| = ∣∣∣∣(Ip(t) − Ipb) ∫ t t−τ G′(s)ds ∣∣∣∣ = ∣∣∣∣ ∫ t t−τ [−Sg(G(s) −Gb) −SiG(s)(Ii(s) − Iib) −SiIib(G(s) −Gb)](Ip(t) − Ipb)ds ∣∣∣∣ = ∣∣∣∣ ∫ t t−τ [(Sg + SiIib)(G(s) −Gb)(Ip(t) − Ipb) + SiG(s)(Ii(s) − Iib)(Ip(t) − Ipb)]ds ∣∣∣∣ 28 J. JIN, J. LI, R. XU, L. YU, AND Z. JIN ≤ 1 2 (Sg + SiIib) (∫ t t−τ (G(s) −Gb)2ds + τ(Ip(t) − Ipb)2 ) + 1 2 SiGM (∫ t t−τ (Ii(s) − Iib)2ds + τ(Ip(t) − Ipb)2 ) = 1 2 [ R ∫ t t−τ (G(s) −Gb)2ds + K ∫ t t−τ (Ii(s) − Iib)2ds + τ(R + K)(Ip(t) − Ipb)2 ] . (A.4) It follows from (A.3) and (A.4), we derive that dV1(t) dt ≤ −A1R(G(t) −Gb)2 −A1 K 2 [ (G(t) −Gb)2 + (Ii(t) − Iib)2 ] | −B1(E + dpi)(Ip(t) − Ipb)2 + (B1 + C1)E(Ip(t) − Ipb)(Ii(t) − Iib) −C1(E + dii)(Ii(t) − Iib)2 + B1 L 2 [ (Ip(t) − Ipb)2 + (G(t) −Gb)2 ] + B1 LR 2 ∫ t t−τ (G(s) −Gb)2ds + B1 LK 2 ∫ t t−τ (Ii(s) − Iib)2ds + B1τ L(R + K) 2 (Ip(t) − Ipb)2 = − ( A1R + A1 K 2 −B1 L 2 ) (G(t) −Gb)2 −B1 [ E + dpi − τ(C + D) − L 2 ] (Ip(t) − Ipb)2 − [ C1(E + dii) + A1 K 2 ] (Ii(t) − Iib)2 + (B1 + C1)E(Ip(t) − Ipb)(Ii(t) − Iib) + B1C ∫ t t−τ (G(s) −Gb)2ds + B1D ∫ t t−τ (Ii(s) − Iib)2ds, (A.5) where C = LR 2 ,D = LK 2 . Second, for V2(t), it follows dV dt ≤ − ( A1R + A1 K 2 −B1 L 2 ) (G(t) −Gb)2 −B1 [ E + dpi − τ(C + D) − L 2 ] (Ip(t) − Ipb)2 − [ C1(E + dii) + A1 K 2 ] (Ii(t) − Iib)2 + (B1 + C1)E(Ip(t) − Ipb)(Ii(t) − Iib) + B1Cτ(G(t) −Gb)2 + B1Dτ(Ii(t) − Iib)2 = − ( A1R + A1 K 2 −B1 L 2 −B1Cτ ) (G(t) −Gb)2 −B1 [ E + dpi − τ(C + D) − L 2 ] (Ip(t) − Ipb)2 − [ C1(E + dii) + A1 K 2 −B1Dτ ] (Ii(t) − Iib)2 + (B1 + C1)E(Ip(t) − Ipb)(Ii(t) − Iib). A NOVEL IVGTT MODEL INCLUDING INTERSTITIAL INSULIN 29 Clearly, if the conditions (H1) and (H2) hold, then dV dt < 0. Therefore, the steady state of Eb is globally asymptotically stable according to the Lyapunov Theorem. This completes the proof. � Jun Jin, Complex Systems Research Center, Shanxi University, Taiyuan, Shanxi, 030006, P. R. China. Jiaxu Li (corresponding author), Department of Mathematics, University of Louisville, ON, USA. Email address: jiaxu.li@louisville.edu Rui Xu, Complex Systems Research Center, Shanxi University, Taiyuan, Shanxi, 030006, P. R. China. Lei Yu, Shanxi Key Laboratory of Mathematical Techniques and Big Data Analysis on Disease Control and Prevention, Shanxi University, Taiyuan, Shanxi 030006, China. Zhen Jin (corresponding author), Complex Systems Research Center, Shanxi University, Taiyuan, Shanxi, 030006, P. R. China; Shanxi Key Laboratory of Mathematical Techniques and Big Data Analysis on Disease Control and Prevention, Shanxi University, Taiyuan, Shanxi 030006, China. Email address: jinzhn@263.net 1. Introduction 2. Formulation of the novel IVGTT model adapted from a physiological model 3. Global stability of the equilibrium of the novel model 4. Studies of the time courses X(t) and Ii(t) with IVGTT data 4.1. Data 4.2. Method of generating the time courses 4.3. Results 4.4. Physiological meaningful parameter values of E and would not destabilize the equilibrium 5. Discussions Acknowledgements References Appendix A. Details of the proof of Theorem 3.2