A comparison among three maximal mathematical models of the glucose-insulin system

The most well-known and widely used mathematical representations of the physiology of a diabetic individual are the Sorensen and Hovorka models as well as the UVAPadova Simulator. While the Hovorka model and the UVAPadova Simulator only describe the glucose metabolism of a subject with type 1 diabetes, the Sorensen model was formulated to simulate the behaviour of both normal and diabetic individuals. The UVAPadova model is the most known model, accepted by the FDA, with a high level of complexity. The Hovorka model is the simplest of the three models, well documented and used primarily for the development of control algorithms. The Sorensen model is the most complete, even though some modifications were required both to the model equations (adding useful compartments for modelling subcutaneous insulin delivery) and to the parameter values. In the present work several simulated experiments, such as IVGTTs and OGTTs, were used as tools to compare the three formulations in order to establish to what extent increasing complexity translates into richer and more correct physiological behaviour. All the equations and parameters used for carrying out the simulations are provided.


Introduction
Type 1 diabetes mellitus is a pathological condition in which blood glucose levels reach excessively high values because of absent or insufficient insulin production [1][2][3]. It occurs mainly in juveniles, due to a combination of genetic determinants and environmental factors [1][2][3]. The worldwide spread of the disease has reached the 8.5% of all diabetes cases in the world, with about 210,000 affected children in the US only [4,5]. Type 2 diabetes mellitus (T2DM) individuals, instead, exhibits excessive blood glucose concentrations mainly due to the insufficient efficacy of circulating insulin to stimulate tissue glucose uptake (insulin resistance) and is correlated with obesity [6,7]. T2DM affects more than 400 million people around the world and it is estimated to double within the next 10 years [8,9], with the increasing in obesity, age and the population size of ethnic groups at higher risk [6]. While metformin, troglitazone, sulfonilureas and other oral hypoglycemic drugs are, together with an improved lifestyle with increased a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 exercise and weight loss, the first treatment choice for most patients with Type 2 diabetes [8], individuals with Type 1 diabetes mellitus require exogenous insulin administrations to supply the lack of endogenous pancreatic insulin production. Insulin cannot be administered orally and must be delivered subcutaneously, either with discrete boli (multiple daily injections, MDI) or through insulin pumps. Recently, much research efforts has been addressed to the development of the so-called "artificial pancreas" (AP) [10][11][12], a system that includes an insulin pump connected to a continuous glucose monitor (CGM) and to a (possibly model-based) control algorithm: the pump delivers an insulin amount calculated on the basis of predicted glycemia. In this case, insulin is administered via a small cannula inserted through the skin into the subcutaneous tissue. Insulin is released as a basal continuous infusion at pre-programmed variable rates, which differ from individual to individual, and as insulin boli at mealtimes. Different strategies, using different mathematical models needed to describe and predict the dynamics of glucose and insulin, exist for the development of the control algorithms embedded in the APs [13][14][15][16]. The accuracy of these models is important in successfully determining the optimal therapy, i.e. the appropriate insulin dose to deliver in order to avoid both hyperglycaemic and hypoglycaemic episodes. The algorithms to be embedded in the pumps must therefore be developed, in-silico, testing their performance over a wide variety of possible physiological situations against a good mathematical model of the glucose-insulin system, implemented in software. The availability of physiologically plausible mathematical models to be used as part of the control algorithm and as a tool for simulating and testing the system in silico is therefore critical for the development of safe and effective insulin infusion control algorithms. It is mandatory that these models are able to reproduce the correct physiological behaviour under different conditions and/or experimental procedures. It is to be noted that many mathematical models of the glucose-insulin system exist. Some of these may be relatively simple approximations, to be used in order to interpret specific clinical testing procedures and identifying specific parameters of interest from relatively small data sets. Others may be relatively complex representations of what is known about the different aspects of glycemic control, including meals, insulin and other hormones, several distribution compartments and so on. We will henceforth denote the simpler models as "compact" and the more extended, hypothetically complete models as "maximal". What is needed for control algorithm development is therefore a physiologically correct maximal model.
In the present work we analyse and compare three such maximal models of the glucose/ insulin system, all three of which have been used in the development of infusion control algorithms: the Hovorka model [17][18][19], the Sorensen model [20,21] and the more recent UVAPadova model [22][23][24][25][26][27]. The UVAPadova implementation used in the present work refers to the S2017 model version [24], which is the latest updated version of the model.
The implemented version of the Sorensen model derives from the equations presented in his original PhD thesis [20], with the corrections introduced in [28] and modified to account for both (possible) subcutaneous administration of insulin and for the reduced insulin-sensitivity, which would be expected in a diabetic patient with respect to a normal individual [1]. The three models present different levels of complexity and a comparison among them gives the opportunity to understand to which extent increasing complexity translates into richer and more correct physiological behaviour.

Comparison among models simulations
In the present work three maximal models of the glucose/insulin system are compared based on glucose and insulin predictions following several simulated experiments. The models analysed are the Hovorka model [17][18][19], the Sorensen model [20,21] and the UVAPadova model [22][23][24][25][26][27]. The Sorensen model is the most complex; it has been extensively described in the Author's PhD thesis and has been thoroughly analyzed in a previous work [28], correcting some errors which were present in the original description by the Author. An extended version of the model, including a representation of the gastrointestinal tract, was also presented in [28], and is used in the present work to simulate experiments where glucose is orally administered. Sorensen's model consists of three sub-models (one for glucose, one for insulin and one for glucagon) that describe the time-course of the variable concentrations in the brain, liver, heart and lungs, periphery (tissue and muscles), gut and kidney. In its original version, it includes pancreatic release of insulin, which in the present work is not considered, given our primary attention to applications for Type 1 Diabetes Mellitus (T1DM). For completeness, all equations are reported in the Appendix. The set of parameter values adopted by Sorensen in his original work is compatible with the normal physiological response to any type of simulated perturbation experiment. The peculiar conditions of T1DM are approximated by Sorensen introducing modifications to the original model. This modified version is the one adopted in the next section. On the contrary, the Hovorka model and the UVAPadova model were originally formulated to represent the physiological behaviour of T1DM individuals, and no modification is necessary for the purpose of the present work. The comparison among the three models was conducted in three steps: • a series of in-silico experiments were set-up and the three models were compared in terms of insulin and glucose concentrations over time; • the three models were adapted to observations of glucose and insulin from an Oral Glucose Tolerance Test performed on a normal individual: the procedure allowed the estimation of the amount of insulin that would be needed to be administered as a bolus in order to obtain the observed time courses, together with some other crucial model parameters. The estimation procedure followed a weighted least squared approach, with weights ω i (i = 1, . . ., n) the inverse of the squared expectations. The variance-covariance matrix of the estimates was obtained with a linear approximation of the model at the optimum by computinĝ s 2 ðĴ TŜÀ 1Ĵ Þ À 1 , where J is the Jacobian; σ 2 × S is the variance-covariance matrix of the observed vector; S is a diagonal matrix with elements (i, i) equal to the squared expectations; s 2 is calculated as 1 nÀ p P iôi � ðy obs i Àŷ i Þ 2 , with n the number of observations and p the number of free parameters. The symbol^is used to indicate quantities computed at the optimum.
• the estimates obtained in the previous step were used to simulate subsequent two other OGTTs, with the administration of the same amount of glucose and insulin given in the first OGTT, with the aim of simulating glucose and insulin concentrations during one day as to mimic three meals.
The in-silico experiments. The implemented in-silico experiments are described as follows: • Intra Venous Glucose Tolerance Test (IVGTT): a continuous administration of basal insulin of 6.67 mU/min, in conjunction with 0.5 g/kg of glucose administered over 3 minutes.
• IVGTT + insulin bolus: a continuous administration of basal insulin of 6.67 mU/min, in conjunction with 0.5 g/kg of glucose administered over 3 minutes, accompanied by a bolus of 1000mU delivered in 1 minute.
• OGTT + basal insulin administration: an oral administration of 100 g of glucose in 1 minute was simulated together with a continuous administration of basal insulin of 6.67 mU/min.
• OGTT + basal insulin + insulin bolus: an OGTT of 100 g of glucose was administered over 1 minute with a continuous basal insulin delivery of 6.67 mU/min in combination with 1000 mU IVITT bolus in 1 minute.
All the above experiments were performed by setting the initial conditions (which also represent the steady state conditions) at a glucose concentration of 5 mmol/L with a basal insulin delivery of 6.67 mU/min for each model considered.

The Sorensen T1DM model
The original version of the Sorensen model was adapted by the Author himself to represent a subject with T1DM. In this process of adapting normal physiology to impaired physiology, the model was modified by removing the pancreatic insulin secretion sub-model and fixing the scale of absolute concentrations of the metabolic source and sink functions in such a way that the diabetic response to any combination of circulating glucose, insulin and glucagon concentrations would have been the same as that of normal individuals subjected to similar conditions. The model adapted in this way therefore represents what would be a "normal" response in a subject with T1DM: for this reason it suffers from an important limitation in that it does not take into account the physiological abnormalities typically present in association with diabetes. In fact, the Author stated that the objective of modelling diabetes condition in the context of his work was ". . .to provide a basis for designing and assessing improved insulin therapies, and in particular for developing an insulin infusion algorithm for closed-loop insulin delivery based on blood glucose measurement.". In this perspective, the comparison of the efficacy of different therapeutic regimens might be considered as largely independent of the details of the physiological model adopted: this is of course no longer true when comparing different models with one another.
A set of normal glucose and insulin concentrations at baseline (Table 1) must be adopted to calculate metabolic rates during diabetic model simulations. This stems from the fact that the post-absorption steady state in the insulin-treated diabetic subjects cannot be determined by the model parameters themselves (as is the case for simulations of normal subjects) as it is dependent on an external forcing function (γ IVI ), which represents the input rate of peripheral venous insulin administration (the therapeutic regimen).
In the case of T1DM modelling, an iterative method for variable initialization must be adopted. Assuming that the response of diabetic subjects to circulating glucose concentrations is the same as in the normal situation, and setting the steady-state values of local glucose and insulin concentrations to values compatible with normal physiology, the glucose concentrations that are at equilibrium with the basal insulin imposed by the external variable γ IVI must be calculated. The procedure is performed by Sorensen and is reported in the Appendix in the subsection "T1DM Sorensen Model initialization".
The equations derived from the steady state conditions are the same as those obtained under normal conditions (model for normal subjects) except for the equation related to I H0 , where the r PIR function is set to 0 and the γ IVI is included in the numerator: I J0 ¼ I H0 ð6Þ All the equations of the Sorensen model are reported in the Appendix (subsection "The Sorensen Model"). As mentioned above, the γ IVI input was introduced by Sorensen into the equation for I H (insulin heart and lung compartment). This obviously represents a simplification, because in the treatment of patients with T1DM insulin is administered subcutaneously, via bolus or continuous infusion. To make the three models comparable, two subcutaneous compartments were therefore added to the revised Sorensen model (which already includes the gastro-intestinal compartment [28]). The model of the subcutaneous compartments used is equivalent to that present in the Hovorka model. The formulation adopted in the UVAPadova model is marginally more complicated: if the parameter k a1 is set to zero and the parameter k a2 is set to same value as the parameter k d , then the two formulations are equivalent. All the parameter values and descriptions are reported in Table 2.

The Hovorka model
The Hovorka model includes two equations for glucose kinetics (amount of glucose in plasma and tissue). Input into the plasma compartment is determined by endogenous glucose production, which depends on plasma insulin concentration, and by absorption through the gastrointestinal compartment. The equations and the values of the model parameters are shown in the Appendix. The term related to Endogenous Glucose Production (EGP), which appears in the final part of Eq [165], represents a linear inverse relationship between glucose production and insulin. This formulation could lead to negative EGP predictions in the presence of high levels of insulin concentrations or in correspondence with particular values of some parameters: a representation that incorporates a saturated effect would be more realistic. The gastrointestinal tract is represented by two compartments: the absorption glucose compartment    (D1), which is fed by the glucose equivalents of ingested carbohydrates, and the conversion compartment (D2) through which uptake of glucose occurs through a linear transfer to the plasma compartment. Insulin is released into the body by means of two subcutaneous compartments, S 1 and S 2 , and then into the bloodstream which represents the plasma insulin compartment. All the values of parameters are shown in Table 3.

The UVAPadova model
The UVAPadova model used in the present work is the S2017 formulation presented in [24], which is an updated version of the original 2013 model [23]. This new formulation includes two new routes of insulin administration: inhaled insulin and intradermal insulin. Some parameters (k p3 , V mx , k p1 ), which in the 2013 version were assumed to be constant, in the later version are made to be time-varying functions. Since their formulations have not been reported in the original work, in the present study we assumed piecewise constant functions for k p3 and V mx from the inspection of Fig 2 in [24], whereas k p1 was set to a constant value. In addition, a new variable (k ir ) is also added, which represents a decreasing factor of insulin dependent glucose utilization (U id ). These latest modifications were included into the model to account for variability of metabolism over 24 hours. The UVAPadova S2017 version makes use of two compartments for glucose (amounts of glucose respectively in plasma and tissues) and two compartments for insulin (amounts of insulin in plasma and liver). The glucose enters the system from the liver (EGP) and the gastro-intestinal tract [22]. Glucose exits the system due to renal elimination and due to glucose utilization, which in turn is divided into two terms: the U id function, insulin-dependent utilization (whose correct formulation is reported in the 2013 version [23]) and the constant U ii , uptake of glucose by the brain and erythrocytes. Insulin appears in plasma via three routes of administration: subcutaneous, inhaled and intra-dermal. The subcutaneous insulin sub-model is described in [25], while the intra-dermal model appears in [29] and the inhaled model is presented in [26]. The action of insulin on glucose is delayed by the introduction of variables that act both on EGP (decreasing as insulin increases) and on insulin-dependent glucose utilization U id (increasing as insulin increases). Once again it should be noted that at high values of insulin concentration, the UVAPadova model would predict negative EGP's. The glucagon sub-model is composed of only one differential equation including: • endogenous glucagon production (input); • subcutaneous glucagon administration (input); • glucagon elimination (output).
Glucagon affects glucose by enhancing its production (EGP increases with increasing glucagon levels in plasma); this effect also occurs through a delayed action. All the values of parameters and their sources are shown in Table 4.

Results
A first set of simulations showed that without making any change to the Sorensen model parameter values, the model is insufficient to predict reasonable time courses of blood glucose concentrations in diabetic subjects undergoing an OGTT: the curve reaches a maximum glycemia of 10 mM, a much lower value than that observed with the other two formulations (18 mM) and presumably observed in such patients; also, the time required for plasma glucose to return to its basal value is approximately 300 min, about half the time required for the Hovorka and UVAPadova models. One reason for the observed divergences from the expected behaviour relies on the assumption that, apart from the defect in insulin production, the physiology of a diabetic individual is otherwise the same as that of a normal individual. This assumption represents an understandable simplification in the absence of further information, but diabetic people are actually known to suffer from reduced insulin sensitivity as well [1], and avoiding to consider diabetic insulin resistance could lead to misleading results. In Fig 1 it can be seen that glucose concentrations forecasts by the Sorensen model differ substantially from those by the other two models. This suggests the need to modify the values of those Sorensen model parameters involved in the description of insulin-dependent glucose uptake. In order to check whether Sorensen's model was qualitatively different from the other two models, parameter values were estimated for Sorensen's model by adapting its predictions to the corresponding Hovorka and UVAPadova predicted time-courses for an OGTT plus insulin experiment (OGTT + basal insulin + insulin bolus in-silico experiment). In order to obtain comparable predictions, three steps were followed: 1. The peripheral venous insulin volume (V I PV ) of the Sorensen model was increased to initialize steady-state insulinemia to a value as close as possible to the initial insulin concentrations observed for the Hovorka and UVAPadova models under the same basal insulin infusion conditions. Parameters F PIC and F LIC , which represent fractional peripheral and hepatic insulin clearance, respectively, were also increased.
2. Plasma glucose concentrations, after oral administration of 100g of glucose without insulin delivery, were simulated with the Hovorka and Sorensen models and then compared to determine the values of the Sorensen parameters involved in the description of insulinindependent glucose production and elimination. The parameters and functions involved in the two processes are r B HGP (the constant basal glucose production), r JGU (the constant rate of intestinal glucose utilization) and M G HGP (the hepatic glucose production) depending on parameters β 0HGP , β 1HGP and β 2HGP , which were all reduced to make the liver less sensitive to circulating glucose concentrations.
3. In the final step the plasma glucose concentrations from the Sorensen model, following an oral administration of 100g of glucose in combination with both an insulin infusion and an  insulin bolus, were made as close as possible to the Hovorka time courses, by modifying the parameters involved in insulin-dependent glucose uptake (M I PGU ). The modified parameters were β 0PGU and β 1PGU , both decreased in order to reflect reduced peripheral tissue sensitivity to insulin, as expected in diabetic subjects.   Table 5 ("(before)" and "(after)" columns respectively). In this figure the predictions from the Sorensen model are comparable with those obtained under the Hovorka and UVAPadova

PLOS ONE
Three maximal mathematical models of the glucose-insulin system model, exhibiting similar maximum concentrations despite a faster return to baseline conditions. Similar time courses were obtained by slightly modifying parameters F PIC and F LIC but by tripling V I PV that was set to a value very close to that assumed by UVA/Padova. Large changes were also necessary for the parameters involved in the insulin-dependent glucose utilization which were decreased by about ten times. This seems however to be reasonable for a diabetic individual.
The other planned simulations, reported in the subsection "The in-silico experiments", are shown in Figs 3-5. The figures show a similar behaviour of the three models with slight divergences: the Hovorka time courses differ from the trends observed for the Sorensen and UVA-Padova models in the OGTT experiment without insulin administration (Fig 3); Sorensen's predictions deviate from the other two in the IVGTT experiment (Fig 4). In the latter figure, the Sorensen model shows reduced insulin sensitivity compared to the other two formulations, exhibiting a much slower return to the basal values, more in line with the profile of a diabetic individual.

OGTT model fitting
The three models were compared in terms of their ability to adapt to observed glucose concentrations from a normal individual undergoing an OGTT with the administration of 100 g glucose. Data were taken from Sorensen's PhD thesis [20]. The fitting procedure was performed by minimizing the sum of the weighted squared residuals (weighted least-squares estimation, WLS, with weights the inverse of the squared expectations). The choice to use real data from a

PLOS ONE
Three maximal mathematical models of the glucose-insulin system normal individual derives from the unavailability in the literature of OGTT data from diabetic subjects, since OGTT is not a standard procedure performed on diabetic people. It should be underscored that model parameters needed to be assessed numerically via fitting, because the parameter values used in the simulations above (representing the response of a diabetic individual), were inadequate to represent the physiological behaviour of normal subjects, who show different rates of absorption and production of glucose.
For each of the three models the fitting procedure allowed the estimation, among other things, of the amount of insulin administered as a bolus; the basal insulin was instead determined in such a way that all three models started (i.e. at time zero, before the glucose and insulin bolus administrations) from the same level of glucose concentration.
The list of the estimated parameters for the three models, together with their before and after estimation process values, are reported in Table 6. The last two columns of the table report the Standard Deviations (SDs) and the Coefficients of Variation (CVs) of the estimated parameters. Since the obtained estimate of the parameter k b3 was essentially zero, the parameter was set to 0 and no variability for it was computed. CVs larger than 100% are not reported and they are identified only as being >100%. This happened for all the free parameters of the UVAPadova model. Fig 6 shows the performance of the three formulations.
For the Sorensen and Hovorka models, the parameters left free to vary are those related to the external insulin input (γ SCIin1 and u 1 , respectively), to the insulin sensitivity mechanism (some parameters in the M I PGU function for the Sorensen model and parameters k b1 and k b3 for the Hovorka model) and to the transfer rates that appear in the subcutaneous insulin compartments (leaving parameter τ S in Eq [177] to vary only for the Sorensen formulation).
The UVAPadova parameters involved in the fitting procedure are those that appear in the external insulin input (u scBolo ), in the representation of insulin sensitivity (k p3 , V mxB and K m0 ), in the delayed effect of the insulin (k i and p 2u ), in the mechanisms of glucose transport between the plasma and tissue compartments (k 1 and k 2 ) and in the effect of glucagon on glucose production (ξ and k H ).
While for the Sorensen and Hovorka model it was necessary to optimize the values of four and three parameters respectively, for the UVAPadova model it was necessary to leave ten parameters free to vary to obtain a good fit of the model predictions to data. The greater number of parameters to be estimated for the UVA/Padova formulation may be due to its great complexity. While it is true that the Sorensen model includes the largest number of equations, it should also to be noted that all parameters are set to values compatible with normal

Fig 6. A)Blood glucose concentrations of Sorensen (dashed line), Hovorka (solid line), UVAPadova (dashed dotted line) models with data points from Sorensen PhD thesis (asterisks); B) Hovorka (solid line) vs data points (asterisks); C) Sorensen (dashed line) vs data points (asterisks); D) UVAPadova (dashed dotted line) vs data points (asterisks).
https://doi.org/10.1371/journal.pone.0257789.g006 physiology. The parameters of the Hovorka and UVA/Padova models are instead indicated for a diabetic individual, but the more compact formulation of the Hovorka model requires fewer modifications to obtain a good fit.
The estimated boluses of insulin required by the three model formulations are: 56387.7 pmol (8886.6 SD, 15.76% CV) (approximately 8 IU), 2361 mU (402.19 SD, 17.03% CV) (approximately 2.3 IU) and 849 pmol/kg (2.508e+05, >100% CV) (approximately 8 IU) for the Sorensen, Hovorka and UVAPadova model respectively. While a similar amount of insulin is necessary for the Sorensen and UVAPadova formulations, the Hovorka model estimated a value four times lower than that estimated by the other two. Fig 6 shows the glycemia time-course predicted by the three models after the fitting process. As expected, the Sorensen model produces the best fit of predictions to observed concentrations. The Hovorka and UVAPadova models seem to be insufficient to represent the final part of the experiment: the Sorensen model is able to predict the rebound observed around the minute 200, where after a decrease below the basal conditions, there is a recovery towards the baseline.
Subsequent three OGTTs in one day. Fig 7 shows the results obtained with three subsequent OGTTs in one day (at 7:00, 12:00 and 19:00) with parameter values set to the estimates obtained by adapting the models to the OGTT data in Fig 6. The Sorensen model appears to be the only one capable of reproducing exactly the same glycemic pattern in the three sub-experiments, with a return to pre-bolus conditions after each OGTT. The UVAPadova formulation is the one for which subsequent OGTTs bring glucose concentrations to lower and lower values (until they fall below a 2mM glycaemia). This feature could be overcome by adopting decreasing values for the parameters V mxB and k p3 , expressing peripheral and central insulin sensitivities respectively, as observed in normal individuals [30]. However, Hinshaw [30] demonstrated that there was no evidence of differences between breakfast and dinner in terms of glucose disappearance in Type 1 diabetic subjects. In the present simulated experiment the values of the two aforesaid parameters were kept constant during the day at the values obtained in the fitting procedure. The Hovorka model predicts glucose concentrations lower than those observed in the first sub-experiment both in the second and third OGTT, nevertheless never producing concentrations below 3mM.

Discussion
Much work has been done within the scientific community, and is still being done, on the study of appropriate models of the glucose/insulin system, aimed at supporting the development of algorithms for controlled and automatic administration of insulin (the "artificial pancreas"). These models must be able to correctly describe the relevant physiology and need to be identified on each single individual: the ability of a model to provide reliable predictions of the glucose and insulin time courses allows the development of robust control algorithms for automatic glucose control in the management of T1DM patients. Among the models present in the literature, the Sorensen model [20], the Hovorka model [17][18][19] and the UVAPadova model [22,[24][25][26][27] are most frequently used to represent virtual patients to this end. The Sorensen model appears to be, among these three, the most complete and detailed in terms of physiological description and parameter values, with its 22 nonlinear differential equations and 135 parameters. Conversely, the UVAPadova model is the most recent and, judging from the number of publications citing it, is the most frequently used, also because its 2013 [23] version was approved by the FDA.
One limit in using this model, however, lies in the difficulty in deciding the values of its many parameters (about 100). In fact, although the mathematical description of the model appears complete from the aggregated publications describing it, the values of several of its parameters are not published, and this prevents the use of the model by potential users. From a usability point of view, the Hovorka model is the simplest, with few model parameters, and therefore easier use for simulation purposes. With the aim of making them available to the interested scientific community, this work provides a complete description of all three models, together with the values of all of the respective parameters.
The present work compares the three models in terms of their performance when simulating the response of a T1DM individual to different glucose stimuli under different types of insulin administration. The comparison among the three models was performed by implementing four types of in-silico experiments. • The Sorensen [20] model can be defined as a maximal model, as it describes the relationship and interactions between the most important actors of glucose metabolism (glucose, insulin and glucagon) for both a normal individual and a diabetic. This means that, compared to the other two, it provides a description of the role of the pancreas, so as to allow not only the simulation of a Type 1 diabetic patient but also of a Type 2 diabetic subject, as well as a normal individual, for which insulin secretion is still maintained or fully guaranteed, respectively. As mentioned above, the model is well documented with regards to the description of both equations and parameters. It does not provide for a mathematical formalization of the gastrointestinal tract and subcutaneous compartments of insulin, but these two drawbacks can be easily overcome by recovering the missing parts from other model formulations.
• The UVAPadova 2018 [24] model is a maximal model that describes the relationship among glucose, insulin and glucagon only for patients with Type 1 diabetes. Compared to the Sorensen model, the UVAPadova model includes a description of gastric emptying and glucose absorption and offers the possibility of considering different types of insulin and glucagon administration routes (subcutaneous, intradermal, inhaled). Most of the model equations are well documented, but some of these are described only in qualitative terms (as for the k p1 and k ir functions in Eq [119] and in Eq [123] respectively): these time-varying parameters were therefore kept constant throughout the simulations. The most important flaw however lies in the lack of values for some parameters. Table 4 reports the values of the parameters used in the in-silico experiments, column Ref; some of them have been found in the literature and the sources are provided (listed with numbers), some have been determined (indicated with the letter D), some have been set at known or reasonable values (denoted by the letter F), the rest of the unknown parameters have been calibrated (and are indicated with the letter C). Therefore, the model, on the basis of what is reported in literature, is not of immediate implementation and use. The present work provides the scientific community with the most complete description having appeared so far of this model's equations and parameters, gathering all the information available from the several sources in the literature [22,[24][25][26][27].
• The Hovorka model [18] is the simplest among the three analysed models; it describes the relationship between glucose and insulin in subjects with in Type 1 diabetes: therefore, like the UVAPadova model, it does not present a description of the secretion and the release of endogenous insulin. While simpler than the other two, this model does provide a clear formalization of gastro-intestinal absorption after oral administration of glucose, and it includes subcutaneous compartments for the representation of external insulin administration. The Hovorka model appears to be easy to understand, sraightforward to implement, and well documented [17][18][19].
The Sorensen model, when used for a diabetic individual, is insufficient to adequately describe the response to any type of experiment if no adjustment is made in terms of parameter values. This is due to the fact that, since people with Type 1 diabetes are at risk of severe hyperglycaemia when undergoing perturbation experiments, no perturbation data are available from the literature, so for this Author it was not possible to derive parameter values under these altered conditions. Sorensen therefore adopts the same formalization and quantification used in normal individuals to approximately describe the physiological behaviour of a diabetic subject, apart from the exclusion of the representation of insulin secretion. Fig 1 highlights the behaviour of the Sorensen model, clearly different compared with the other two. The Hovorka and UVAPadova models instead, show similar time courses, reaching the same maximum value with a slight difference in the return to basal conditions. After making the appropriate modifications to the Sorensen model, as described in the "Methods" section, the resulting predictions resemble those obtained with the other two formulations (see Fig 2), although a faster recovery towards the basal values is observed. It is however clear that, in the absence of actual observation data on diabetic subjects, it is not possible to state which of the three models comes closer to describing the correct physiology. Fig 3, mimicking an OGTT without insulin bolus, shows similar trends for Sorensen and UVAPadova, with similar maximum values reached. The Hovorka model predicts higher glucose concentrations and a slower return to baseline conditions. The higher concentrations could be due both to higher than expected endogenous glucose production and to lower tissue glucose uptake in the absence of an insulin bolus.
In Fig 4, where the IVGTT experiment is simulated in conjunction with an insulin bolus, the Hovorka and UVAPadova models show a similar time course, while the predictions of the Sorensen model appear to be qualitatively different. Although no observations on actual patients are available to demonstrate the plausibility of the three models, the predictions obtained with the Sorensen model seem to be more in line with the expected response of a diabetic subject since exhibit a slower return towards pre-experiment values.
In the IVGTT experiment without insulin bolus administration (Fig 4) UVAPadova and Hovorka are quite in agreement, while Sorensen shows a glucose trend more consistent with an insulin-resistant profile. These results have been obtained by modifying some model parameters of the Sorensen formulation, in particular those relating to the description of both central and peripheral insulin sensitivity (Table 5). It is likely that,by changing the values of other parameters, a greater similarity of the predictions of the three models may be obtained. However, in the absence of clear evidence pointing to the better performance of one model with respect to the others, introducing changes in the parameter values chosen by the Author after a thorough literature search is not advisable unless supported by physiological justifications.
The adaptation of the three models to real OGTT data of a normal individual emphasizes the ability of all three models to adapt to real glycemic trends. Even if the setting in which the comparison is made is not optimal (observations are made on a normal individual and not on a diabetic patient), this procedure allows us to investigate the ability of the models to adapt to real data, leaving some parameters free to vary and estimating the administered insulin dose necessary to reproduce the glucose observations. A "good" model should be able to elucidate the relationship between glucose and insulin and should be able to predict a recovery to baseline conditions, with a time course as close as possible to that of a normal subject, with a reasonable amount of insulin. This aspect is important when these models are used in model-based control algorithms, so that the ability of a model to adapt to real observations becomes an essential feature and deviations from what is observed emphasize important physiological deficiencies. Fig 6 shows a very good performance of all three models in adapting to the data, but we can see that the Sorensen model is the only one able to predict the recovery phase with a rebound after a hypo-glycaemic period produced by the administration of insulin, reaching the last available data point. This could be due to the fact that the parameter values used for the Sorensen model were those originally adapted to the physiology of a normal individual. Surprisingly, while the Hovorka model tries to predict a recovery phase, which will indeed occur at a later time, the UVA/Padova model appears to stabilize at a lower blood glucose levels despite a greater number of free model parameters (ten for UVa-Padova, four for Sorensen and three for Hovorka). The estimated amount of insulin administered as a bolus, needed to obtain the predictions of Fig 6 are 2.3 IU, 8.5 IU and 8 IU for the Hovorka, UVAPadova and Sorensen formulations respectively. The much lower value observed for the Hovorka model could depend on the estimated values of the other free parameters: insulin sensitivity, explicitly represented in the Hovorka model by the parameter k b1 , is estimated in fact at about 1.6 × 10 −4 , and if insulin sensitivity is high, then necessary amount of insulin to be delivered decreases. In this context however it should be noted that a lower k b3 (roughly representing central peripheral insulin) determines greater production of endogenous glucose, which on one hand compensates for the augmented sensitivity to peripheral insulin, while on the other causes final recovery to attain a higher blood glucose concentration. Concurrently with a large bolus of insulin obtained for the Sorensen model, the parameters of the M I PGU function (representing insulin sensitivity) all decreased by about a factor of two. For the UVAPadova model an external insulin input (u scBolo ) of approximately 8 IU is required. A much lower K m0 parameter (approximately zero) in Eq (123) is consistent with much improved insulin sensitivity, accelerating the effect of insulin on glucose utilization.

Conclusion
The three models seem to reproduce all the simulated experimental situations quite well without obvious divergences. The Sorensen model produces predictions similar to those of the other two models once some parameter values are modified, suggesting that, with suitable adaptation, this model could be used to also represent the physiology of diabetic subjects. An updated formulation including both the gastrointestinal tract and the subcutaneous insulin deposit compartment should be used in this case.
With the information available at the present time, the final choice about which model to use lies in the confidence that the experimenter places on how plausibly the mathematics represents the underlying physiology, as well as on the simplicity, robustness and versatility of the formulation. A more complex model with a large number of parameters might in principle fit better with observations, but the complexity of a model not only makes identification statistically harder, but it suffers from possible over-fitting and consequently fragile, unreliable forecasts. From this point of view, the Sorensen model is the one with the greatest number of parameters (adding the fixed, determined and free ones). However, since they are well documented and an in-depth bibliographic research has been carried out by the Author, most of the parameters are set to known values, so as to require the estimation of a lower number of free parameters than the UVAPadova model, which in any case is poorly documented. Conversely, Hovorka's model is the simplest, and still fits sufficiently well the observed data. We notice however that the Sorensen model is the only one capable of predicting the rebound phase in the OGTT experiments, where the other two models fail.
From purely in silico experiments it is not possible to draw definitive conclusions on which model is physiologically more credible. The next logical step in the evaluation of these and possibly other maximal models of the glucose-insulin system would be to compare their predictions against actual observational data, obtained with different experimental set-ups in patients with a range of normal and diabetes conditions.

The Sorensen Model
Mass Balance-Glucose BRAIN HEART AND LUNGS Metabolic Source and Sinks-Glucose Mass Balance-Glucagon Metabolic Source and Sinks-Glucagon Glucagon V Γ = 11310ml 1: Sorensen PhD thesis [20] Determined parameters.
a ¼ 5 2½Dose 0 ð1 À bÞ� ð152Þ b ¼ 5 2ðDose 0 cÞ ð153Þ Renal excretion of glucose Insulin effect on distribution/transport of glucose Insulin effect on glucose disposal Insulin effect on EGP released from liver Oral CHO intake expressed as glucose equivalent With M wG ¼ g mol � � glucose molecular weight Glucose in the absorption compartment Glucose in the conversion compartment The glucose absorption rate Amount of short-acting insulin Insulin concentration Insulin absorption rate into the blood Determined parameters.