Consistency of compact and extended models of glucose-insulin homeostasis: The role of variable pancreatic reserve

Published compact and extended models of the glucose-insulin physiologic control system are compared, in order to understand why a specific functional form of the compact model proved to be necessary for a satisfactory representation of acute perturbation experiments such as the Intra Venous Glucose Tolerance Test (IVGTT). A spectrum of IVGTT’s of virtual subjects ranging from normal to IFG to IGT to frank T2DM were simulated using an extended model incorporating the population-of-controllers paradigm originally hypothesized by Grodsky, and proven to be able to capture a wide array of experimental results from heterogeneous perturbation procedures. The simulated IVGTT’s were then fitted with the Single-Delay Model (SDM), a compact model with only six free parameters, previously shown to be very effective in delivering precise estimates of insulin sensitivity and secretion during an IVGTT. Comparison of the generating, extended-model parameter values with the obtained compact model estimates shows that the functional form of the nonlinear insulin-secretion term, empirically found to be necessary for the compact model to satisfactorily fit clinical observations, captures the pancreatic reserve level of the simulated virtual patients. This result supports the validity of the compact model as a meaningful analysis tool for the clinical assessment of insulin sensitivity.


Introduction
Mathematical models for representing and simulating insulin/glucose metabolism, both in normal conditions and under perturbation tests, have been developed and used since the 1960's [1]. Their importance resides in the possibility to better understand the pathophysiology and development of Type 2 Diabetes Mellitus (T2DM) and of its pre-conditions, such as Impaired Glucose Tolerance (IGT) or Impaired Fasting Glucose (IFG). More extended models are also used to simulate virtual patients in order to test control algorithms, (i.e. [2][3][4][5][6] and references within) in situations where experimental procedures are very expensive or cannot be a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 performed for ethical or practical reasons. Models can be more or less complex according to how many sub-components of the entire metabolic system are simulated (glucagon, adrenalin, Free Fatty Acids, etc.) and according to the level at which the relevant physiology is modelled (e.g. whole-body or cellular). Comprehensive models of the glucose-insulin system [7,8] allow to simulate healthy subjects and T2DM patients for clinical research purposes, e.g. for testing control algorithms in silico.
Two test procedures in general use for the estimation of insulin sensitivity are the Intra-Venous Glucose Tolerance Test (IVGTT), which needs to be analyzed by means of a suitable mathematical model [9][10][11]; and the Euglycemic Hyperinsulinemic Clamp (EHC) [12], considered the 'gold standard' for the assessment of insulin resistance, which yields a measure of insulin sensitivity by direct averaging of the final glucose infusion rate. The standard IVGTT is simpler to perform than the EHC, has no significant associated risks and, if appropriately modelled, provides important information about the negative feedback regulation of glucose and insulin in a specific subject. Efforts have been made therefore to develop better models for the interpretation of the IVGTT experimental data set, exhibiting behaviour compatible with physiology (for example bounded and positive solutions of the model, as for example glucose predictions that never go to infinity or assume negative values) [13] and providing stable (assuming similar magnitude among different subjects), precise estimates of the structural model parameters [10]. Minimal models (minimal in the number of equations and in the number of parameters to be estimated) have therefore been developed to estimate the insulin sensitivity of a specific subject from a relatively non-invasive test procedure such the IVGTT. The so-called "Minimal Model" (MM) [9], which is still the most widely used compact model in the clinical setting, was demonstrated [11,13] to suffer from a number of drawbacks, among which poor parameter estimation, with very large parameter coefficients of variation, translating into overall non-identifiability of the model parameters, including in particular the S I index of insulin sensitivity. Moreover, the estimation procedure typically used, decoupling the feedback and estimating separately the two glucose-insulin control arms, provides misleading results as is discussed in greater detail elsewhere [10]. An alternative model of the compact model class, previously published by the present authors [10] and referred to as the Single Delay Model (SDM), was demonstrated to have mathematically consistent solutions, admitting the fasting state as its single equilibrium point and converging back to it from the perturbed state. This model was simultaneously fitted to glucose and insulin IVGTT observations on a heterogeneous population composed of lean, overweight, obese and morbidily obese subjects and was proven to outperform the MM over the whole range of conditions considered [11]. An extension of the SDM was subsequently proposed [14], where in order to model the fate of per os glucose, a gastrointestinal tract model was added to the plasma insulin and glucose dynamic equations and the glucose rate of appearance was derived by describing the absorption of glucose along a sequence of three gut compartments.
With the aim of interpreting a wide range of heterogeneous experimental results related to pancreatic insulin secretion, an islet population model was proposed in 2010 [15], incorporating the original idea of Grodsky of a population of independent controllers, coupled only by circulating glucose levels [16]. This model explains the effect of the islet response to varying glucose concentrations by means of a simple second-order nonlinear model, of the same functional form for all islets, with a random distribution of parameter values over the large number of islets considered. The population of pancreatic β-cells, collected into Langerhans islets, can be viewed therefore as a set of independent, similar, but not identical controllers (firing units) with distributed functional parameters. The islet equations are coupled to a metabolism submodel to complete the description of the feed-back control of the glucose/insulin system dynamics. This islet population model was shown to reproduce very closely a wide array of actually observed, diverse in vivo and in vitro experiments, including the pioneering work of Grodsky, with the same set of working parameters [17]. While the model does not include any dependency on the rate of change of glycemia, it is able to reproduce accurately the double phase of insulin release during a prolonged glucose stimulus: a first phase of impulsive insulin release, immediately upon glucose administration, and a second phase of more gradual release, also linked with the potentiation effect of persistent hyperglycemia on the secretory units.
During the course of the development of the SDM in 2007 it was appreciated that, while the inclusion of some model components was irrelevant to the quality of data fitting (e.g. the inclusion of "glucose effectiveness" or an esplicit representation of the delay of insulin action on glucose), other elements, such as the delay τ and the exponent γ, were indispensable, and neglecting them led to the inability of the model to adequately fit the observations. Parameter τ represents the delay with which the pancreas changes secondary insulin release in response to varying glycemia, while the exponent γ represents the rapidity with which the insulin secretion rate reaches its maximum value with increasing glycemia. It was proven in [10] that both the delay τ and the exponent γ included into a nonlinear sigmoidal function were indeed necessary to interpret data. Now, a few years after the publication of the SDM model, a comprehensive model for insulin secretion exists [15,17], which is in fact able to explain mechanistically why a host of morphologically diverse insulin secretion responses to glycemic perturbations occur. We can therefore go back and attempt to understand why the nonlinearity γ was so essential in the representation of insulin secretion within the compact SDM model: this is the aim of the present work.
The reason why the islet population model is used as a comparison model is that it mirrors the actual anatomic structure of the endocrine pancreas and, at the same time, is able to convincingly reproduce observed insulin secretion patterns after many different types of glycemic perturbations.
Starting from the extended and validated model [17], the following procedure is implemented: a) a sample of IVGTT data sets is generated from the extended model portraying a group of virtual patients with a range of insulin sensitivities; b) the generated IVGTTs are fitted with the SDM compact model; c) a comparison is made between the set of generating extended model parameter values and the set of estimated compact-model parameter values to clarify their relationships.

Compact Model
The delay model used in the present work, henceforth Compact Model (see Fig 1), has been shown to be able to represent well the glucose and insulin concentrations observed during an Intra-Venous Glucose Tolerance Test (IVGTT) [10]. Using this model it is possible to estimate the insulin sensitivity of a patient by fitting the patient's IVGTT data, and this estimation is precise, reproducible and robust [11]. This model has been shown to admit mathematically consistent solutions with physiologically plausible parameter values [18].
Let G(t) and I(t) be, respectively, the plasma glucose concentration [mM] and the serum insulin concentration [pM], while their distribution volumes are indicated, respectively, as V G [L/kgBW] and V I [L/kgBW]. The differential equation representing variation of plasma glucose concentration is: with initial conditions: The glucose basal concentration G b [mM] is the glycemia level before the bolus injection, while G Δ [mM] represents glycemia increase following the bolus. D g [mmol/kgBW] is the intravenous dose of glucose administered at time 0 during an IVGTT experiment. The equation representing the dynamics of plasma insulin concentration I [pM] is the following: with initial condition: The novelty this model introduced is in the second term of (3), representing second-phase insulin delivery from the β-cells. Its functional form is consistent with the hypothesis that insulin production is limited, reaching a maximal rate of release T max ig =V I by way of a Michaelis-Menten or a sigmoidal Hill dynamics according to whether the γ value is 1 or greater than 1 respectively. The constant values T gh and T max ig can be directly computed from the steady state conditions of Eqs (1) and (3) respectively. See the original references for more details [10,11]. Glucose Plasma concentrations arise following an IntraVenous Glucose Dose. Input to plasma Glucose compartment comes also from the liver which produces glucose by means of gluconeogenesis and glycogenolysis processes, while elimination occurs by an independent (first order elimination rate) and dependent-insulin mechanism (second order tissue uptake). High plasma glucose concentrations enhance insulin release in plasma, while insulin is cleared with a first order elimination rate.

Extended Model
The mathematical model of the glucose-insulin system presented in [17], henceforth Extended Model, has been the first model of whole-body insulin secretion to provide a unified explanation of an array of diverse clinical experimental procedures. In fact, the model is able to reproduce, with the same set of (meta)parameters: the low-frequency ultradian oscillations appreciable in the insulinemic signal when a constant enteral feeding is administered to a patient [19]; the entrainment of insulinemia to glycemia when a patient undergoes different Intra-Venous (I.V.) glucose administration patterns with different frequencies [20]; high-frequency insulinemia oscillations triggered by I.V. administration of very small amounts of glucose [21]; reproduction of the glycemia and two-stage-insulinemia (first and second phase) curves, upon simulated I.V. glucose administration during an IVGTT experiment. In this latter case, moreover, the model is able to mimic the clinically observable glycemia and insulinemia curves both for Normal Glucose Regulation (NGR) and for different pre-morbid and morbid conditions, such as Impaired Fasting Glycemia (IFG), Impaired Glucose Tolerance (IGT), IFG+IGT and Type 2 Diabetes Mellitus (T2DM).
While this model is extensively discussed elsewhere [17], we briefly summarize its features in the following.
The basic paradigm of this model (see Fig 2) is that in the pancreas a multitude of similar, but not identical, independent controllers react to the sensed plasma glucose, which acts as the single "coupling" signal. While the qualitative behavior of the firing units is the same, each unit reacts differently to glycemia, and these heterogeneous performances yield the characteristic insulin responses to different stimuli. Note that the physiological identification of the firing unit could be the β-cells scattered in the pancreatic Langerhans islet, or, by choosing a different level of model granularity, subcellular granules, or, conversely, collections of synchronized βcells within the islets of Langerhans (for more details, see [17]). The pancreas secretory units (circles) release at different times their packet of insulin J n (t) (depending, for the n-th controller, on the threshold B n (t) and on the potentiation level D n (t)) in such a way that at any given time a total quantity J(t) (sum of the amounts J n (t)) of insulin flows into the portal vein (Q 1 compartment), then to the liver (Q i+1 . . .Q L compartments) and finally reaching the plasma (compartment I), stimulating the uptake of glucose by tissues. Glycemia (G compartment), which is raised by glucose hepatic production k 3 , stimulates the production of insulin, closing the loop. A firing unit (unit n) releases at time t � its stored packet of insulin J n (t � ) [pmol/kgBW] whenever circulating glycemia exceeds the current threshold value B n (t � ) [mM], thereupon entering a (relative) refractory state, where further stimulation fails to elicit the release of new hormone (the time-course of the threshold function B n (t) is in general different for each unit n); the refractory state of unit n is represented in the model by instantaneously increasing the glycemia threshold of that unit to a high level R n [mM] whence, over time, it exponentially decreases towards its resting threshold value G n [mM], G n < R n . As soon as the controller fires, its threshold abruptly increases its value and the controller enters a refractory state. From this moment onwards, the threshold value decreases exponentially. When the threshold B n reaches again values comparable with the current glycemia, the unit n is ready to release a new packet again. The size of the packet of insulin itself depends on prevailing glycemias: in fact, the well-known phenomenon of potentiation occurs, that is the ability of the pancreas to respond with progressively increasing insulin amounts to identical glucose stimuli, when these are repeated in close proximity over time [16,22,23]. The variation in size of the packet of insulin (subject to potentiation) for each firing unit is described by the equation for D n (t) [pmol/kgBW] [17].
Each different firing unit n is then characterized by the same triple of equations (B n (t), D n (t), J n (t)), with different parameter values. As described in depth elsewhere [17], their values are randomly sampled from given (usually lognormal) distributions.
The differential equation associated to B n is: where α n [min −1 ] is the rate of recovery of sensitivity of the secretory unit, G(t) [mM] is the external glycemia sensed by all the secretory units (that is the glucose plasma concentration), and δ(�) is a Dirac delta term specifying instantaneous increase of the threshold to the refractory level R n , associated with discharge of insulin, at any time the glucose stimulus G(t) exceeds the controller threshold B n ; χ is the characteristic function of its argument set. The parameter G n is the only one that is not extracted from a lognormal distribution, but instead from a distribution with the following density (see for more details [17]): depending on the two parameters g 1/2 [mM] and ν [#]. The Insulin Secretion Rate (ISR) is defined as: where N is the total number of firing units. Basically, Eq (7) states that the ISR at time t is given by the sum of the insulin packets fired by the units whose threshold B n (t) is smaller than or equal to the current glycemia G(t). Insulin then flows to the portal vein and the liver according to: where Q [pmol/kgBW] refers to the insulin amount in the portal vein/liver compartment. Then, insulin mass Q enters the plasma insulin distribution space according to the following equation: where I(t) [pM] is the serum insulin concentration. Finally, glucose plasma concentration is described by: where: 1. the first term in (10) describes approximately the (supra-threshold) driving glycemia for urinary glucose elimination:ũ 2. the second term is the insulin-dependent glucose uptake, with k 2 [/min/pM] being the rate of glucose uptake by tissues per pM of serum insulin concentration; 3. the third term k 3 (t) [mmol/kgBW/min] refers to the net balance between hepatic glucose output and insulin-independent zero-order glucose tissue uptake (essentially by the brain), with V G [L/kgBW] the apparent distribution volume for glucose.
In order to account for the noisy time-course of the last term k 3 , a stochastic model has been employed: where � k 3 [mmol/kgBW/min] is a central value for k 3 (t) and ξ(t) [mmol/kgBW/min] is a stochastic process. For further details, see [17].

Modifications of the two models
In the present work, we use the fact that the Compact Model [10] is able to fit IVGTT data of a virtual patient generated by the extended model [17]. In a previous work [24] we obtained very good fits of the Compact Model on the Extended Model using directly parameters reported in [17] for generating virtual patients with profiles ranging from NGR to T2DM; nevertheless, we introduce here some modifications of the Extended Model in order to better reproduce observed IVGTTs from a larger sample of clinical tests.
In particular, since both the Extended and the Compact Model do not present any term of insulin-independent glucose elimination, we modified the glycemia Eqs (1) and (10) respectively as: and Consequently from (13), the new value for T gh is computed at the equilibrium as: It is to be noted that the Compact Model does not include the glycosuria term, which is instead present in the Extended Model-see Eq (10).

Virtual patient fitting
Before producing virtual patients with the Extended Model, some fine tuning of the model parameters was thought to be advisable in order to better match glycemia and insulinemia trends of Normal Glucose Regulation (NGR) patients undergoing an Intra-Venous Glucose Tolerance Test (IVGTT), as reported in the literature (original data which served to set the best parameters values derived from Panunzi et al. [10], from which only the subsample of lean subjects was considered). A few modifications in the parameter/metaparameter values of the Extended Model have therefore been carried out, as reported in Table 1: as it may be noticed, the glucose net balance mean value ( � k 3 ) has been halved, while the mean value of the rate of recovery of sensitivity of the secretory unit (α n ) has been decreased, obtaining a slower insulin dynamics. Moreover, the renal glucose elimination threshold G u has been increased according to the literature. In Fig 3 are reported as thin black lines the 5th and 95th percentile of the data analyzed in [10], while the median of 100 virtual NGR patients generated by the Extended Model with the new parameter values reported in Table 1 is shown as a thick red line. This trend is coherent with NGR glycemia and insulinemia curves observed during IVGTTs as reported in the literature [25][26][27].
With this new set of (meta-)parameters, virtual patients exhibiting different degrees of disease (ranging from a Normal Glucose Regulation (NGR) patient to a Type-2 Diabetes Mellitus (T2DM) patient) were generated by modifying parameters as reported in Table 2. In the same way as previously done [17], we obtained virtual T2DM patients by lowering peripheral insulin-sensitivity (parameter k 2 ), increasing liver glucose production � k 3 (reflecting lower hepatic insulin sensitivity), by lowering the size of the insulin packets (parameters μ(ρ n ) and mð � D n Þ), and by lowering firing thresholds (whose distribution depends on the parameter g 1/2 ). The resulting pancreatic behavior still shows a higher insulin production given by the higher number of recruited firing units (due to lower firing thresholds), in the effort to compensate the lower insulin sensitivity and under the hypothesis that these T2DM patients are still able to compensate. However, the ability of these subjects to further increase their pancreatic insulin secretion responding to higher glycemias is severely limited. We generated IVGTT glycemia and insulinemia curves for a total of 500 virtual patients, divided into 25 degrees of morbidity, ranging from "pure" NGR to "pure" T2DM condition in 25 uniform steps (see Table 3); this choice is a good tradeoff between the computational time required to generate all the cohort and the value necessary to have an "average" patient for   [27][28][29] are comparable to the ones generated by the Extended Model. Observation data were then extracted for each virtual patient from the glycemia and insulinemia curves at times −10, −5, −0.5, 2, 4, 6, 8, 10, 15, 20, 25, 30, 40, 50, 60, 90, 120 min. The "observations" were finally fitted with the Compact Model, estimating the following vector of parameters θ: Since the simulated glucose bolus has been set in such a way to produce a sudden glycemia increment of 10 mM for every patient, G Δ in Eq (2) was fixed to 10 mM as well. The glucose distribution volume V G was set to 0.2 L/kgBW and kept fixed during the generation of all virtual patients.
The fit for each virtual patient was performed by solving a constrained optimization problem (lower and upper bounds were considered) using a Nelder-Mead algorithm minimizing the following loss function λ: where G i and I i are observed glycemia and insulinemia at time t = i, whereasĜ i andÎ i are the corresponding estimated glycemia and insulinemia. This function was adopted in order to consider at the same time both glycemia and insulinemia data, using a normalization factor.

Results
Virtual NGR patients were generated by using the same model parameter values employed for producing curves reported in Fig 3. Fig 4 reports the result of the fitting procedure on data extracted from a virtual NGR patient, while Figs 5 and 6 report the fitting performance for a virtual pre-diabetic patient and for a virtual T2DM patient respectively; Table 2 reports the Compact Model estimated parameters compared with the generating Extended Model parameter values.
Since the insulin-sensitivity parameter is present in both models (k XGI for the Compact Model and k 2 for the Extended Model), we reported in the same plot (Fig 7)  We recall that the relationship between circulating glycemia and insulin secretion for each virtual patient depends on the threshold distribution of the parameters G n (see Eq (5)), randomly extracted from the distribution f(g) in Eq (6). The shape of this function depends on the pathophysiological condition of the virtual patient and therefore on the value of its g 1/2 parameter (see Table 3).
We now define z 2 [0, 1] [#] to be a coefficient of pancreatic reserve, representing, for each given subject, the fraction of total firing units able to react to an abrupt raising in glycemia above the baseline. It can be computed as: where G b [mM] is the basal glycemia of the patient (before the injection of the glucose bolus), ν = 2.5137 [#] (see [17]) is constant for all patients and g 1/2 [mM] varies, passing from NGR to T2DM, according to Table 3.  healthy (k XGI � 1 � 10 À 4 , blue square), pre-diabetic (0:5 � 10 À 4 �k XGI < 1 � 10 À 4 , red triangle) and diabetic (k XGI < 0:5 � 10 À 4 , black dot). In this figure, healthy and pre-diabetic patients cluster into two tight distinct sets; diabetic patients are scattered over a different, widely diffused set. The three sets present a clear differentiation based on the value of the pancreatic reserve z: healthy subjects exhibit a high pancreatic reserve ranging from 0.7 to 0.9, pre-diabetic patients are characterized by a lower reserve value ranging from 0.4 to 0.8 and diabetic subjects present with a very low index of pancreatic reserve, below 0.4. NGR and pre-diabetic patients exhibit estimated values of γ which lie in the neighborhood of a falling parabola on the plane z-γ, reported in Fig 11 as a blue curve, their γ increasing as the pancreatic reserve z decreases. The diabetic patients, instead, although having low values of z, exhibit widely different estimated values for γ. The behavior of the z-γ relationship for normal and prediabetic subjects is clear. In the case of diabetic patients, the actual pancreatic reserve is so low, and the maximal possible increase in insulin secretion is so small, that the γ parameter from the Compact Model is difficult to identify precisely, with the corresponding worsening of the scatter of the estimates as z goes to lower and lower values. Fig 12 shows, for the 500 virtual patients, the relationship between the estimated γ and the value of T max ig for the Compact Model, representing the maximum reachable Insulin Secretion Consistency of compact and extended models of glucose-insulin homeostasis Rate. The value T max ig has been computed as follows: whereK xi [min −1 ] is the estimated constant disappearance rate for insulin, I b is the basal insulinemia of the patient, V I = 0.25 [L�kgBW −1 ] is the constant distribution volume for the virtual patient (from the Extended Model) and G � = 9 [mM] is the glycemia at which the insulin secretion rate is half of its maximum and is fixed for all patients. Both V I and G � values are the same as in [10]. Also for Fig 11, when plotting T max ig against γ, the subjects cluster again in three well defined groups: NGR patients (blue squares) exhibiting low values of γ, present a T max ig which is typically below 4 pmol/min/kgBW; pre-diabetic patients (red triangles) exhibit higher values of γ and higher T max ig ; diabetic patients (black dots), showing values of γ in a wide range (due to the above-mentioned difficulty of precisely estimating γ when increments of insulin secretion are limited), present a T max ig above 3 [pmol/min/kgBW]. This pattern is consistent with the modest needs for increased insulin secretion in healthy, insulin-sensitive subjects, in whom, on the other hand, no ceiling for insulin secretion is readily apparent, hence the low nonlinearity coefficient γ, and vice versa for insulin-resistant subjects. Fig 13 reports the dependency of T max ig on z for each patient. We may observe the same clear separation among the three different classes of patients apparent in the other plots. Consistency of compact and extended models of glucose-insulin homeostasis It can be seen that NGR patients exhibit a high value of pancreatic reserve z and low values of T max ig , since they do not need high insulin production due to their good insulin sensitivity; conversely, pre-diabetic and diabetic patients, with a lower pancreatic reserve, have a higher T max ig , since they need to compensate their lower insulin-sensitivity.

Discussion
In the past 50-60 years many attempts at mathematical modelling of different aspects of the glucose-insulin system have been made. Some contributions concentrated more on the insulin secretion mechanisms [15,16,30], others more with short-term modelling of perturbation experiments [9,10,13,14,31], and some publications describe more complex models, from multi-organ models [8,32] to maximal models for in-silico Type 1 diabetes virtual patient simulation [33][34][35][36][37]. Our group has been active in this area for many years. We have addressed in particular total body or organ-level integrated mechanisms, modelling both short-term experiments [10,11,13,14,18], pancreatic insulin release mechanisms based on firing units [15,17,24], long term compensation and disease development [38], incretin effect [39]. Our overall goal is to eventually obtain a coherent suite of models addressing different relevant clinical issues (e.g. IVGTT and OGTT for insulin sensitivity determination, or dynamics of insulin secretion), suitable for integration into an overall scheme. These models should be individually effective in dealing with the specific dynamics they approximate, should exhibit appropriate qualitative behavior from a mathematical viewpoint, should have desirable statistical properties leading to a-posteriori parameter identifiability and should also be consistent with one another, shedding light on the actual mechanisms involved in glucose homeostasis.
In particular, the clinical use of a compact, identifiable model should be justified by our confidence that it does capture the basic relevant physiology, if not all possible features of interest. In the present work we are trying to establish the relationship between an extended, validated model (which by its very size cannot be fruitfully employed for the analysis of limited experimental data e.g. from an IVGTT) and a Compact Model of proven effectiveness in estimating insulin sensitivity from small datasets. By generating virtual patients from the Extended Model and fitting them with the Compact Model we are able to link variations in the (simulated) physiology with diagnostic parameter value changes.
Some modifications of the Extended Model parameter values (with respect to previously published values [17]) were found to be necessary in order for simulated IVGTT curves to lie Consistency of compact and extended models of glucose-insulin homeostasis squarely within 90% confidence bounds obtained from sets of clinically recorded IVGTT's from NGR subjects. These minor changes in the parameter values, however, did not affect the reproducibility of the diverse patterns generated from the whole array of experimental procedures, such as the slow (ultradian) [19] and fast oscillations [21], and the entrainment phenomenon due to sinusoidal glucose infusions [40], which were all re-checked with the modified parameter set.
Moreover, in order to align the two models, an insulin-independent glucose elimination term was introduced in the glycemia equation of the Extended Model. The insulin secretion formalization remains, of course, radically different between the complex Extended and the simple Compact Models.
Data of 19 patients from previous clinical studies [11] were used to build 90% confidence envelopes of the observed glycemia and insulinemia time courses. Fig 3 shows that for the set of chosen parameters the median of 100 virtual NGR patients undergoing an IVGTT experiment lies within the 90% envelopes from real observations highlighting the ability of the Extended Model to reproduce well real IVGTT data.
Starting from the NGR parameter values, progressive degrees of clinical worsening (up to T2DM) were simulated by changing some model parameters, such as insulin sensitivity (k 2 ), average net hepatic glucose production ( � k 3 ), and other parameters related to the mechanisms of insulin production (average ρ n , average � D n , and g 1/2 ). The Compact Model was then fitted onto the IVGTT data generated with the Extended Model, and some of the Compact Model estimated parameters were compared with the corresponding values of the generating parameters. Consistency of compact and extended models of glucose-insulin homeostasis Figs 4 to 6 clearly show the good performance of the Compact Model to fit the Extended Model curves, highlighting its ability to reproduce both first and second phases of insulin release. It is evident that as healthy conditions are altered to pre-diabetic and then to diabetic status, the second phase of insulin release appears to be more and more prominent, probably in order to compensate a higher insulin resistance translating into a lesser degree of glucose uptake and persistently high plasma glucose concentrations.
The consistence and robustness of the Compact Model is evident from Fig 7 where the relationship between the two model parameters k 2 (for the Extended Model) and k XGI (for the Compact Model), both quantifying insulin sensitivity, approximates very well the identity line in the k 2 -k XGI plane.
The overestimation of the insulin sensitivity index (k 2 ) for very high insulin resistance patients is due to the absence, in the Compact Model, of the term representing glucose renal elimination, which is instead present in the Extended Model. While, in principle, this structural difference could be absorbed by a lower net HGO (T gh ), this compensation mechanism is insufficient (Fig 9) as it is glucose-independent: this causes a higher glucose uptake, which in the model is represented by a second order term (in glycemia and in insulinemia). Figs 8 and 10, compared with the corresponding Figs 7 and 9, show how the introduction of the glycosuria term in the Compact Model restores the linear relationship between k XGI and k 2 even at low insulin sensitivities, as well as solving the underestimation problem of T gh with respect to � k 3 (the net HGO in the Extended Model). Consistency of compact and extended models of glucose-insulin homeostasis Of greater interest is however the observation portrayed in Fig 11. When pancreatic reserve is maintained (z > 0.4, see Eq (18)), there is a clear (inverse) relationship between said pancreatic reserve and the nonlinear behavior of the pancreatic insulin release represented by parameter γ. When pancreatic reserve is exhausted, the increase in insulin secretion is too small with respect to the increase in glycemia with resulting very variable values of the nonlinearity exponent γ. Fig 13 shows that the maximum insulin production T max ig is small in athletes or very healthy subjects (the pancreas in these subjects is not typically stimulated very much since insulin sensitivity is excellent), increases with BMI and prediabetes (consistently with clinical observations of high insulinemias in these subjects) and then decreases as the clinical picture worsens and (relative) pancreatic insufficiency develops; again consistently with clinical observations.

Conclusion
A refinement of our understanding of the population-of-controllers interpretation of the insulin secretion behavior of the pancreas (in the several degrees of impairment between NGR and T2DM) explains why the Compact Model for the IVGTT needs to explicitly incorporate the Consistency of compact and extended models of glucose-insulin homeostasis saturating nonlinearity of insulin secretion with increasing glycemia, and the fact that this nonlinearity worsens with progressive worsening of the clinical picture.
This consistency in the behavior of the two models, developed for and fitted on very different sets of data, increases our confidence that both models, each in its own domain, are reliable representations of the actual physiology.
It is remarkable that the Compact Model, with only two parameters relative to the insulin secretion mechanism, and estimable on the small datasets from IVGTT experiments, is able to capture a series of detailed physiological mechanisms explicated and represented in the Extended Model: the end-result of a very complex set of events, such as the potentiation and synchronization of the whole β-cell population, along with the introduction of a heterogeneous behaviour of the secretory units (from the distribution of thresholds to the distribution of the rate of recovery of sensitivity) appears to be well captured by the very simple but specific nonlinear function of the Compact Model.

Author Contributions
Conceptualization: Andrea De Gaetano.