MATHEMATICAL ANALYSIS OF A MODEL FOR GLUCOSE REGULATION

Diabetes affects millions of Americans, and the correct identification of individuals afflicted with this disease, especially of those in early stages or in progression towards diabetes, remains an active area of research. The minimal model is a simplified mathematical construct for understanding glucose-insulin interactions. Developed by Bergman, Cobelli, and colleagues over three decades ago [7, 8], this system of coupled ordinary differential equations prevails as an important tool for interpreting data collected during an intravenous glucose tolerance test (IVGTT). In this study we present an explicit solution to the minimal model which allows for separating the glucose and insulin dynamics of the minimal model and for identifying patient-specific parameters of glucose trajectories from IVGTT. As illustrated with patient data, our approach seems to have an edge over more complicated methods currently used. Additionally, we also present an application of our method to prediction of the time to baseline recovery and calculation of insulin sensitivity and glucose effectiveness, two quantities regarded as significant in diabetes diagnostics. 2010 Mathematics Subject Classification. Primary: 92B; Secondary: 92C60.


84
FESSEL, GAITHER, BOWER, GAILLARD, OSEI AND REMPALA 1. Introduction.Diabetes affects approximately 26 million U.S. adults and children, including 7 million undiagnosed cases.An additional 79 million individuals have pre-diabetes, with elevated blood glucose levels below diagnostic cut points [17].Pre-diabetes is a state of abnormal glucose homeostasis in individuals without diabetes, where a deficiency or resistance to insulin is apparent [22,28].Individuals with pre-diabetes are of particular interest to researchers and clinicians because progression to diabetes can be prevented or delayed in this population with early identification [11].
Current standard of practice primarily utilizes three biomarkers for the identification of pre-diabetes: (1) plasma fasting glucose, (2) post-load plasma glucose, and (3) glycated hemoglobin (HbA1c) [2].Each of these tests has distinct limitations [3,12,27] and clinical cut points are selected somewhat arbitrarily.Exact onset time of what is termed "diabetes" is not a discrete event and therefore recommended clinical categories are based primarily on the gradient of the association of these biomarkers with prevalent retinopathy [19,20] and evidence from clinical trials demonstrating that lowering these values can reduce microvascular complications [36].The intravenous glucose tolerance test (IVGTT) can provide important information in particular about first-phase insulin response to glucose [9,18] and can help researchers better characterize individuals at increased risk for developing type 2 diabetes [23,26].
The exact point where diabetes onset occurs remains unknown.A normal pancreatic β-cell is highly adaptable to changes in insulin action.For example, as insulin action is decreased, insulin secretion increases to meet the demand imposed by the change in insulin action.When adaption of the β-cell is no longer sufficient for a given degree of insulin insensitivity, impaired glucose tolerance or type 2 diabetes develops.Insulin resistance-where the muscle, adipose tissue, and liver cells do not use insulin properly-develops when the biological effects of insulin are abnormal for glucose disposal in skeletal muscle as well as endogenous glucose production suppression [33].
Multiple mathematical models of glucose-insulin interaction have been proposed (see, e.g., the review articles [1,15]), but perhaps the most widely known and widely used was suggested by Richard Bergman, Claudio Cobelli, and colleagues by the early 1980s.Their minimal model proposed in [7,8] utilizes a system of ordinary differential equations (ODEs) to represent the joint effect of insulin secretion and sensitivity on glucose tolerance [4].The full ODE model is composed of two subsystems: the first describes glucose clearance via the equations relating glucose and interstitial insulin and the second describes plasma insulin action.The minimal model as well as other similar ODE systems are used as investigatory modeling tools to measure insulin secretion and insulin sensitivity [5,6] in individual patients.The typical methods used to parameterize these models (i.e.MINMOD [31], the Matlab routine gluc mm mle.m [34], etc.) rely on numerical ODE solvers, however, and solutions of such do not yield a straightforward dynamical time course for the glucose concentration, which is perhaps a major reason why the ODE based methods for glucose-insulin dynamics are still not very popular in the experimental diabetologists' community [32].
To address this issue, in the current paper we derive an explicit formula for the glucose subsystem in the non-dimensionalized or scale-free Bergman-Cobelli minimal model [29].We do not change the structure of the non-dimensionalized minimal model, but rather we develop a novel method to solve for the glucose time course.This approach allows us to definitively parametrize individual glucose patterns by using the IVGTT data that inform the assessment of insulin sensitivity.Our method also permits the development of results, such as return time to baseline, that are not possible with current ODE minimal model solvers [14,34].By explicitly solving for the glucose-time relationship for a given individual, we have developed a model which is not only easier to fit with patient-specific data but also significantly more transparent regarding the time evolution of the glucose concentration, a benefit likely to be appreciated by diabetes practitioners.
The paper is organized as follows.In the next section (Methods), we outline the derivation of our explicit formula for the glucose subsystem, first recalling some basic facts about the minimal model and its scale-free extension.In the subsequent section (Results), we show how the derived glucose function may be fit to the IVGTT data of individual patients, and we provide some potential applications of the model.We also compare our model fit with that of MINMOD and two additional more recent models of glucose-insulin dynamics of De Gaetano and Arino [21] and Li et al. [25].We conclude (Discussion) by reviewing observations from our analysis of the minimal model and by suggesting future investigations using our framework.

Model analysis.
In this subsection we present our mathematical formulas for glucose regulation based on solving the minimal model.We begin by recalling some essentials of the Bergman-Cobelli minimal model [7,8].Extending the ideas of Nittala et al. [29], we then shown how one may non-dimensionalize the minimal model in order to simplify the relative glucose-insulin interdependencies.
We then present our explicit solution for the glucose concentration G(τ ) and the interstitial insulin X(τ ).This is accomplished using an intermediary function f (τ ) to link G and X, and we also obtain an implicit formula for the plasma insulin concentration I(τ ).Lastly, we use the initial and long-term conditions for G and X to further characterize the intermediary function.
2.1.1.Model review and non-dimensionalization.Denoting dimensional variables with asterisks, let G * be the glucose concentration, I * be the concentration of plasma insulin, and X * be the insulin in the interstitial fluid [5].The equations governing glucose disappearance during an IVGTT are given by the minimal model [7,8]; namely, where insulin is taken to be a known input function.G b and I b are the baseline glucose and insulin concentrations, respectively, and p * i are various rate parameters.We assume the initial conditions for this system to be: G * (0) = G 0 , X * (0) = 0, and I * (0) = I 0 .Additionally, the long-term values for our variables are: G * (∞) = G b , X * (∞) = 0, and I * (∞) = I b ; that is, the glucose and insulin concentrations return to baseline after a sufficient amount of time has passed.Bergman, Cobelli, and colleagues also developed an equation to describe the IVGTT insulin kinetics given a known glucose forcing function.These two subsystems, one for glucose clearance and the other for insulin action, were not meant to serve as a closed, coupled system for G, X, and I [8], and we will not treat them as such here.
Moreover, due to the difficulties in insulin response detection, we only focus on the subsystem controlling the glucose clearance in this study.
We non-dimensionalize equations (1) in accordance with Nittala et al. (choice c) [29] so that the scaled variables are Assuming max(G * ) = G 0 and max(I * ) = I 0 , both G and I range approximately between [0, 1]. 1 We represent time as τ = t * /t 0 where t 0 = 60 min, and the remaining parameters are then given non-dimensionally as Using these scalings the minimal model becomes where G 1 is a dimensionless quantity: The initial conditions for this scaled system are: G(0) = 1, X(0) = 0, and I(0) = 1, and all three quantities return to zero as τ → ∞.
2.1.2.Explicit solution of the non-dimensional model.As pointed out by Nittala and colleagues [29], X(τ ) can be written in terms of I(τ ) directly.The explicit solution for (3b) is τ 0 e p2s I(s) ds.
However, we may also solve the glucose differential equation (3a) with analytic techniques since X(τ ) is now a known function of I(τ ).
Making the variable transformation Ĝ(τ ) = G(τ ) + G 1 and utilizing the method of integrating factors, we arrive at where From equation (6) we see that X(τ ) can alternatively be expressed in terms of the intermediary function f (τ ); specifically, Note that equation (5) implies in particular that G(τ ) only depends on the insulin level X(τ ) though the intermediary function f (τ ).Further, equating (4) and (7) gives an implicit equation for the relationship between the plasma insulin concentration and the intermediary function: Returning to the proposed initial conditions, we find that G(0) = 1 is satisfied for any choice of f (τ ), but requiring X(0) = 0 dictates that Further analysis on the system, such as applying the conditions as τ → ∞, hinges upon further specifying the form of the intermediary function, f (τ ).
2.2.Implementation.Obtaining an explicit solution for G and X leads to many powerful analytic possibilities.We use our achieved solution to aid in patientspecific model fitting in this paper; however, other applications are conceivable.In what follows we present f (τ ) as a series of exponentials and develop various specifications on this function.We examine a special case of f (τ ) that leads to exponential decay of the glucose and approximate inactivity of the interstitial insulin.We close this subsection by establishing a numerical approach for determining individualized parameters, which depends upon applying our explicit solution to observed human data.

Solution specifications.
In interest of analytic tractability, we assume the intermediary function to be composed of a series of exponentials; that is, let where the parameter set {α i , β i } i=1,...,n ultimately characterizes each patient's glucose-insulin dynamics.The initial condition requirement for X (9) now stipulates that and the long-term condition for X yields max i {β i } = p 1 , since the term with the maximal exponent dominates as τ → ∞.Without loss of generality, we set n i=1 α i = 1.Combining these conditions gives three equivalent quantities: Note that our long-term glucose condition is automatically satisfied as long as max i {β i } > 0, a specification already provided for by the minimal model since p 1 > 0 [7,8].
Next we are tasked with determining the required number of exponential terms to include in the expansion of f (τ ).As outlined in A, it turns out that at least three terms must be included in order to satisfy (11) and to retain the freedom of unique exponents (β 1 = β 2 = . . .= β n ).We therefore choose n = 3 to minimize the complexity of our solution; however, the number of exponentials may be increased if additional conditions are desired.
Again without loss of generality, we pick max i=1,2,3 {β i } = β 3 and use (11) to solve for β 3 .We thus find where Requiring In summary, we find that the minimal model ODE system can be solved analytically where G, X, and I are associated via an intermediary function f (τ ).We write f (τ ) as the sum of three exponentials and determine that the system's initial and long-term conditions require (11).When dictating i α i = 1 and max i {β i } = β 3 , the coefficient of the third term is taken to be while the third exponent is given by (12).Lastly, the remaining coefficients and exponents must satisfy the inequalities (13).Now consider a special case in which all three of our exponents are roughly equivalent, β 1 ≈ β 2 ≈ β 3 .Both expressions in ( 13) hold approximately at equality, and since all exponents are nearly the same, we find f (τ ) ≈ e bτ where b = p 1 ≈ β i for i = 1, 2, 3.The glucose and interstitial insulin solutions (5) and ( 7) reduce to G(τ ) ≈ e −bτ and (15a) Note that if f (τ ) is simplified to a single exponential, then the non-dimensionalized glucose concentration undergoes a simple exponential decay while X(τ ) remains static and approximately zero.The glucose-insulin dynamics in this case are markedly different than those of the full model ( 2.2.2.Computational approach.We now develop a computational approach that utilizes our formulated explicit equations to predict individual glucose dynamics from patient IVGTT data.To this end, we apply an optimization algorithm based on the least-squares error criterion to estimate the coefficients and exponents of the intermediary function along with G 1 .The non-dimensional quantity G 1 can be calculated from G b and G 0 , which are theoretically known for each patient; however, regarding G 1 as a parameter allows for added model flexibility [29].We initialize our method by supplying an individual's IVGTT data as well as starting values for our unknown parameters.Once the glucose data has been nondimensionalized, we employ Matlab's fmincon.m library function for parameter estimation.(In the work presented here, we apply the active-set method, though other schemes could easily be selected.Additionally, we choose the maximum function evaluations to be 2000 and a functional tolerance level of 10 −5 .)Since the fmincon.mroutine allows for both linear and nonlinear constraints, we provide the following conditions: The first requirement (16a) dictates that at least one of the found exponentials be nonnegative.(Recall max i {β i } = p 1 > 0.) The last two stipulations guarantee that β 3 = max i=1,2,3 {β i }.Once parameters that both satisfy the constraints and minimize the glucose error are found, the remaining parameters, α 3 and β 3 , are calculated from ( 14) and ( 12), respectively.Equation ( 5) now completely determines the predicted glucose time course and may be graphed for visual comparison with the empirical data.
After the intermediary function is fixed, we use fmincon.monce more to determine the rate constants p 2 and p 3 by minimizing the difference between the left and right hand sides of the implicit insulin equation (8).The individual's observed IVGTT insulin values serve as I and numerical integration is completed with the trapezoid method.The only necessary constraints for this step are p 2 , p 3 ≥ 0.
For the special case discussed at the end of Section 2.2.1 where β 1 ≈ β 2 ≈ β 3 , we can simply perform a linear fit to the logarithm of our glucose data since G ≈ exp(−bτ ).We use Matlab's polyfit.m routine and have p(τ ) = p1 τ + p2 where p = ln G to determine b = p1 .(Note that allowing for the constant p2 in polyfit.m is analogous to admitting G 1 into our parameter search set.)To avoid log-transforming nonpositive numbers, we set any nondimensionalized glucose measurements below ε = 10 −2 to ε.
We now have natural starting values to run the full three-exponential method; specifically, [α 1 , β This special case sufficiently models the glucose dynamics of some individuals and returns β 1 ≈ β 2 ≈ β 3 .(See, for example, the patient data described in Section 3.2.1.)The simplified fit is not appropriate for many others, however, and new initial values must be provided for {α 1 , α 2 , β 1 , β 2 }.The starting value for G 1 is always taken to be G

Results.
3.1.IVGTT data.We present the results of applying our analytic-numeric formulation to observed human data in this section.For each individual studied, we determine parameters describing the glucose clearance pattern with our minimization routine.We observe that the resolved pattern may follow either a simplified exponential decay or may require our more complex full model.We provide representative examples of each in Sections 3.2.1 and 3.2.2,respectively, and we remark that many other patients likewise fall into one of these two categories.We also compare our results to those obtained using a traditional ODE solver method and two other models of glucose-insulin dynamics based on delay-equations.Lastly, we discuss applications of our approach including predicted time to baseline recovery and the quantification of well-known diabetic composite parameters: insulin sensitivity and glucose effectiveness.We execute our proposed method on IVGTT data gathered from several human patients.Patients included in this research were participants in an intervention study examining the effect of weight loss on high-density lipoprotein functionality in a population of overweight and obese adults with pre-diabetes not on diabetes medications or statins.Data presented in this paper were extracted from the baseline examination.The study was reviewed and approved by the Institutional Review Board at the Ohio State University Wexner Medical Center; all study participants provided written informed consent prior to data collection.
Glucose and insulin concentrations were collected from each individual at 0, 2, 5, 8, 10, 16, 19, 22, 25, 30, 40, 60, 90, 120, 140, 160, and 180 minutes post-bolus injection.The first glucose recording, measured simultaneously with injection, is taken to be the known baseline glucose, G b .The second time point marks the maximal glucose value in all but a few patients in our data set, so our constructed numerical method is used to predict the glucose time course from two minutes onward.

Patient 1:
Single-exponential fit.The glucose clearance of our first example patient can be sufficiently modeled with a single exponential; that is, G(τ ) ≈ e −bτ where b = p 1 ≈ β i for i = 1, 2, 3. Applying a linear fit to the logarithm of our nondimensional data yields a predicted value of b = 1.602.We use this value to initialize the full computational program as specified by (17).The resulting coefficients from this method are listed in Table 1.While the exponential values are not precisely equivalent to b, 2 we do find that β 1 ≈ β 2 ≈ β 3 .
Figure 1 shows the full, three-exponential fit for Patient 1 which is characterized by the parameters given in Table 1.Though not shown in the figure, we remark that the simplified, single-exponential fit closely follows this predicted curve.
3.2.2.Patients 2 and 3: Full model fit.Our second and third example patients exhibit glucose clearance dynamics that cannot be modeled with a single exponential function.Applying the simplified model to the dataset of Patient 2 yields b = 1.908; however, using this value to initiate the full computational program as in (17) gives 2 The minor discrepancy between b and β i is attributed to setting a tolerance for the logarithm values and to the fact that we are using two different error metrics to obtain our predicted parameters.Given N data points, the method for the full model minimizes however, in the single-exponential fit we minimize N k=1 ln G data,k − ln G f it,k 2 .an approximate single-exponential fit that does not match the patient data well, primarily because this glucose pattern returns to baseline from below.Alternatively, we launch our full method with the starting parameters taken as [α 2 ] = [0.5, −10, −1, 2].From these initial values, we arrive at the coefficients and exponents listed in Table 1 for Patient 2. Note that the exponentials found with this method are not equivalent (β 1 = β 2 = β 3 ) and that the full model now appropriately fits the glucose data as represented by the solid red curve in Figure 2 (top).We also predict the interstitial insulin concentration from formula (7), as illustrated by the blue curve in Figure 2 (top).
We must use the full, three-exponential model to achieve an appropriate fit for the glucose patterns of many other patients in our dataset including Patient 3. Following the same procedure described for Patient 2, the final resulting parameters for Patient 3 are listed in Table 1.The model glucose and interstitial insulin curves for Patient 3 are displayed in Figure 2 (bottom).

3.3.
Fit comparison with ODE solver and delay differential equations.We now consider how our analytic method fares when compared against other glucose models, including the minimal model upon which it is based, and two more complex delay-time methods.
Algorithms which fit IVGTT data with the minimal model, including variants of the MINMOD software, traditionally iterate between two steps: estimation of patient-specific parameters, frequently accomplished with the method of least squares, and glucose curve fitting, conventionally performed via a numerics-based ordinary differential equation (ODE) solver [13,31,34].(Recall that while we use least squares for parameter fitting, we do not rely upon an ODE solver because we have developed an explicit time-glucose formula.) An additional refinement for glucose-insulin modeling is to incorporate the baseline insulin concentration as a state-value which pre-exists before the injection.This improvement necessitates delay differential-equation models, of which two of the better-known are those of De Gaetano and Arino ( [21]) and Li et al. ([25]).In this section we compare the fit achieved using our routine to that obtained by solving the minimal model with a typical ODE solver (namely the Matlab routine gluc mm mle2012.m developed by Natal van Riel [34,35]) and also to the fit of De Gaetano's and Li's models (which we apply to our data using the Matlab routine dde23.m).
In Figure 3 we see the curves predicted for Patients 1 and 3's glucose by the four models -our explicit, the ODE minimal-model, and De Gaetano and Li's delay  Glucose model fit (solid red) along with ODE solver gluc mm mle2012.m fit (dashed blue) to the patient glucose data (black dots) for Patients 1 and 3. De Gaetano's and Li's models, fitted to both glucose and insulin data (insulin data not shown), appear in dashed magenta and cyan and coincide so closely as to be indistinguishable.(The fit for Patient 2 is very similar for the three non-ODE models; see Table 2.) differential equations (DDE) models.One obvious initial conclusion is that Li and De Gaetano's models are significantly less accurate than the minimal and exact models in the description of glucose, for these two datasets.(In the case of the second patient, the three non-ODE-models are of very comparable accuracy and hence their graphs are omitted-see Table 2.) We regard this tendency towards inferior fitting as a simple consequence of the two more complicated models' need to account for insulin concentrations.
Another point is that the ODE-solver-based method (i.e. the minimal model, numerically approached) is inferior to the other three models in the initial region t = 0 . . .8.
In Table 2 we quantify the fit of all four models by computing their corresponding coefficients of determination3 (R 2 values) for the three patients' glucose data (which we utilize as our measure of fitness, following De Gaetano's lead [21].Also following De Gaetano, we ignore those glucose measurements recorded before time t = 8.) 4We see that when only glucose dynamics are considered, our explicit model performs better for Patients 1 and 3, and with comparable accuracy for Patient 2.
Finally, we remark upon the deviations between explicit and MINMOD ODE models which occur near the beginning and the end of the IVGTT time course.The main reason for these differences is that van Riel's ODE solver (similarly to the DDE models) uses a weighting scheme which downweights the first few data points; whereas, we have chosen to weight each measurement equally.More on this subject can be found in the Discussion.

Time to baseline recovery.
The time required for a patient's glucose level to return to baseline after bolus injection reflects the individual's health status [10,24,30].Since our approach produces an explicit expression for the temporal glucose concentration, we can easily extend the time course beyond the length of the IVGTT experiment to predict when a patient returned to his/her baseline glucose level.In this study we define the time to baseline recovery (t b ) as the point where the individual's predicted glucose level recovers to within 5% of G b and remains as such for 10 minutes thereafter.
Consider the glucose trajectories of Patients 1 and 2. As seen in the left panel of Figure 4, we predict that Patient 1 returns to glucose baseline 163 minutes after injection, a time point within the duration of the IVGTT.Patient 2, however, does not return to baseline within the confines of the 180-minute test.Our model predicts that Patient 2 returned to baseline at t b = 372 min.This type of simple predictive analysis could possibly shorten the recommended IVGTT duration or serve as a metric to indicate existence or severity of disease.Further note that the ODE solver methods do not yield an explicit time-glucose relationship and therefore cannot be extended to predict time of return to baseline.
3.4.2.Insulin sensitivity and glucose effectiveness.An important measure for health quantification is insulin sensitivity, S I .From the Bergman-Cobelli minimal model [7,8] and the scaling (2), we have that .
Firstly, note that determining the insulin sensitivity from the model fit for Patient 1 is impossible.Since f (τ ) ≈ e bτ and b ≈ p 1 , equation ( 8) is only satisfied if p 3 ≈ 0. Our insulin error minimization technique is consistent with this observation and returns a value for p 3 that is below computational tolerance for Patient 1.This clearly demonstrates a weakness in assuming the glucose concentration follows strict exponential decay.If the glucose is assumed exponential, the minimal model concludes that insulin does not affect the glucose clearance dynamics at all.The glucose effectiveness, on the other hand, can be calculated for all patients during the glucose curve fitting.This diabetic measure, S G , is equivalent to the first rate parameter; that is, S G = p * 1 .We list the resulting dimensional values for {p * 1 , p * 2 , p * 3 } along with the composite parameters {S G , S I } for all three patients in Table 3.We also compare the values obtained with our method (columns labeled "Explicit") to those found using van Riel's method (columns labeled "ODE") [35].The associated units for each parameter are min −1 for p * 1 , p * 2 and S G ; [mU/L] −1 min −2 for p * 3 ; and L•mU −1 •min −1 for S I .

4.
Discussion.The minimal model developed by Bergman, Cobelli, and colleagues nearly three decades ago endures in modern diabetes literature for quantifying insulin sensitivity (S I ), glucose effectiveness (S G ), and for inferring diabetic risk from frequently sampled IVGTT data.In this paper we have provided a mathematical solution to the minimal model thereby translating this coupled differential equation system into explicit, time-dependent formulas which rely on the analysis of glucose dynamics only.Exploiting the straightforwardness of our solution as well as contemporary advances in numerical minimization algorithms, we have resolved some of the issues with the original MINMOD software [31].For example, with our approach it is now possible to administer computational constraints to ensure that p i ≥ 0. We have compared our model fit and parameter estimations to the results of a typical minimal model ODE solver as well as the fits of two more recent models utilizing delay differential equations of De Gaetano and Arino [21] and Li et al. [25].Additionally, we arrived at critical composite parameters (S I and S G ), and we showed how our explicit solution can predict future glucose dynamics including the time of return to baseline (t b ).
In our analysis we identified one basic classification scheme based on whether or not a patient's glucose concentration pattern followed a simple exponential decay.As pointed out in the literature [5], a single exponential curve cannot account for the rich dynamics of glucose clearance which occurs in multiple stages, often with the glucose concentration returning to baseline from below.A single-exponential fit cannot model such scenarios since it only accounts for the intrinsic glucose decay; that is, dG dτ = −p 1 G, which means the insulin action is disregarded and X ≈ 0. Consequently, a single exponential decay pattern may indicate insulin insufficiency just as dips below baseline may indicate insulin overcompensation.Our full, threeexponential model shows promise for capturing a wide variety of clearance patterns, including those suggesting insulin deficiency or surplus, which may aid researchers studying diabetes onset.
In Pacini and Bergman's computational implementation of the minimal model (MINMOD) [31] as well as several other approaches, weighting schemes are commonly adopted, especially for the first few glucose data points which are often omitted since extracellular mixing prevails during this time period.While we did not employ a weighting scheme for the work presented here, we have for comparison attempted 1) relative weighting and 2) downweighting of the first few IVGTT measurements but we found that these refinements did not significantly influence our results.
In our analysis of the IVGTT data the explicit model outperformed decisively the minimal model ODE solver method as well as in two out of three cases also the two more recent alternative models of glucose-insulin control system of De Gaetano and Arino [21] and Li et al. [25].Although these two models are based on more detailed physiological analysis of glucose-insulin interactions and apply more sophisticated modeling tools, when restricted to the glucose subsystem only they do not seem as capable of reproducing the correct IVGTT dynamics as the explicit model.We hope to conduct further, more comprehensive comparisons to this end in the near future.
Indeed, our presented work serves as proof of concept for future larger studies aimed at establishing clearer reference ranges among healthy, pre-diabetic, and diabetic individuals.Apart from potentially simplifying the IVGTT protocol itself, the favorable predictive properties exhibited by the explicit model could be also useful when analyzing datasets with missing or less frequent time point data.With a given mathematical solution to the minimal model, we can potentially reconstruct glucose response curves using a smaller number of observed data points (and no insulin data) as well as exploit analytic techniques to test parameter sensitivities.The ability to predict glucose behavior beyond measured time points affords us an opportunity to better characterize pathophysiological features of the diabetes development process.Additionally, our explicit approach to modeling IVGTT data could conceivably be adapted to the more widely used oral glucose tolerance test (via the model proposed by Caumo et al. [16], for example).and equation (11) becomes Without loss of generality, we assume max i=1,2 {β i } = β 2 .We also have i α i = 1, which for this case means α 2 = 1 − α 1 .Equation ( 18) thus gives and the only way to satisfy ( 19) is to set β 1 = β 2 .This requirement induces the same model inflexibility mentioned for the n = 1 case.We thus find the smallest number of exponentials that can satisfy (11) with multiple unique exponents is n = 3.While three exponentials are so required, we turn our attention to a special case in which all three exponents are approximately equal; that is, β 1 ≈ β 2 ≈ β 3 .Since the sum of the coefficients is one, we find Furthermore, the concentration of insulin in the interstitial compartment (7) remains near zero for all times in this special case: X(τ ) ≈ be bτ e bτ − p 1 = 0.

1 =
G b /(G 0 − G b ), and we choose p .036 • (I 0 − I b )to reflect typical quantities in the insulin error minimization[13].

Figure 1 .
Figure1.Glucose model fit (solid red) along with its 95% Gaussian confidence envelope (dashed red) for the data from Patient 1 (black dots).This individual's predicted glucose clearance roughly follows exponential decay.The dotted, horizontal line represents the baseline glucose concentration and has been included for reference.

1 )Figure 2 .
Figure2.Glucose model fit (solid red) along with its Gaussian 95% confidence envelope (dashed red) to the patient data (black dots) for Patient 2 (top) and Patient 3 (bottom).The glucose clearance depicted here cannot be modeled by a single exponential and requires the use of our full, three-exponential method.The predicted interstitial insulin concentration (solid blue) is calculated from equation(7).The dotted, horizontal line indicates the individual's baseline glucose concentration for reference.

Figure 3 .
Figure3.Glucose model fit (solid red) along with ODE solver gluc mm mle2012.m fit (dashed blue) to the patient glucose data (black dots) for Patients 1 and 3. De Gaetano's and Li's models, fitted to both glucose and insulin data (insulin data not shown), appear in dashed magenta and cyan and coincide so closely as to be indistinguishable.(The fit for Patient 2 is very similar for the three non-ODE models; seeTable 2.)

Figure 4 .
Figure 4. Extended dynamics of the parameterized glucose function.The time of return to baseline glucose (t b ) is indicated by the vertical line; all times thereafter are shaded gray.Our explicit solution predicts t b = 163 min for Patient 1 (left) and t b = 372 min for Patient 2 (right).Each patient's baseline glucose value is indicted by the dashed line for reference.

Table 2 .
R 2 values for goodness-of-fit of four models to glucose data, for times t ≥ 8.

Table 3 .
Dimensional Rate Parameters and CompositesThe glucose concentrations of Patients 2 and 3, however, do not follow simple exponential decay, and while we do not predict the plasma insulin time course with this approach, we can effectively approximate p 2 , p 3 , and thus S I for this group of individuals.