A revised Sorensen model: Simulating glycemic and insulinemic response to oral and intra-venous glucose load

In 1978, Thomas J. Sorensen defended a thesis in chemical engineering at the University of California, Berkeley, where he proposed an extensive model of glucose-insulin control, model which was thereafter widely employed for virtual patient simulation. The original model, and even more so its subsequent implementations by other Authors, presented however a few imprecisions in reporting the correct model equations and parameter values. The goal of the present work is to revise the original Sorensen’s model, to clearly summarize its defining equations, to supplement it with a missing gastrio-intestinal glucose absorption and to make an implementation of the revised model available on-line to the scientific community.


Introduction
Diabetes is a metabolic disorder that affects millions of people worldwide and depends on the reduced capacity or the complete failure of the pancreas to produce insulin, possibly combined with the resistance of tissues to insulin action. People affected by diabetes, if not properly treated, present elevated levels of glucose in the bloodstream (hyperglycemia), and if not well controlled may develop long-term severe chronic complications such as retinopathy, stroke, nephropathy, neuropathy, or even experience dangerous and potentially lethal episodes of hypoglycemia [1][2][3][4][5][6].
In this framework, the need for developing tailored treatment approaches has led to the use of mathematical models able to represent the glucose/insulin dynamic response of the single patient, for instance with the objective of fully automating-insulin delivery. With the goal of maintaining glucose levels inside a narrow range of values, a number of closed-loop control methods have been designed, often relying on a limited or sparse number of observations, approach that is often referred to with the general term artificial pancreas (see e.g. [7][8][9][10][11][12] and references therein); these systems base their efficacy, in terms of in-silico testing, on reliable mathematical models of glucose homeostasis. Over the last sixty years, several efforts have been made for the representation of the glucose/insulin system through the use of compartmental models, some simpler ("minimal models"), with few equations and parameters [13][14][15][16], others more complex, including interactions among many players (insulin, glucagon, free fatty acids‥). Among the most widely used comprehensive multi-compartment models of the glucose-insulin system there are the UVa-Padova model [17][18][19][20], the Hovorka model [21,22] and the model developed by Sorensen [23]. In particular, the Sorensen model is perhaps the most complex among these, incorporating 22 differential equations (mostly nonlinear), representing glucose concentrations in the brain, heart and lungs, liver, gut, kidney and periphery, with about 135 parameters (including the initial conditions of the state variables). The values of the many model parameters were decided on the basis of a careful literature research. Many researchers subsequently based their work on the Sorensen model ( [24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39] represent just an example), but because of its complexity only few authors used it for the development of control algorithms [31,35,38,39]. The model delivers a detailed and well-documented representation of the physiological mechanisms, but, possibly due to its complexity, several scientists who employed and implemented it did not perform a sufficiently thorough, painstaking analysis of the model itself, thereby inheriting some errors from the original Sorensen presentation.
Moreover, while intravenous glucose administration experiments and continuous glucose infusion rate experiments can be explicitly represented by the original model, oral glucose administration lacks such an explicit representation. In fact, the 100g oral glucose test described by the Author in his original work is simulated by means of the direct introduction of a gut glucose absorption rate term in the gut mass balance equation (r oga ), bypassing the glucose pathway from stomach to gut and adopting an empirical functional formulation for the time curve of the absorption rate. In the process of adapting the model to experimental data, the Author adjusted the rate of gut glucose absorption in order to match observed blood glucose concentrations with model-predicted peripheral blood glucose concentrations.
Another limitation of the Sorensen model is the inability of the pancreas sub-model to appropriately secrete insulin in response to an oral glucose load, which would have required some description of the action of incretin hormones. This limitation was bypassed by the Author by empirically estimating the required rate of pancreatic insulin output to adjust predicted insulin concentrations with respect to observed insulin experimental data [23].
The objective of the present work is twofold: first, we want to highlight and correct the errors appearing in the original version of the Sorensen's model, to implement a carefully revised version of the model, and to make the implementation available on-line to the scientific community. The implementation is provided both in user-to-machine and machine-tomachine versions at the address: http://biomatlab.iasi.cnr.it/models/login.php (access as a Guest).
Matlab code is also downloadable from the same link. The second objective is to enrich the model by introducing a description of the gastrointestinal tract (to simulate alimentary glucose intake, digestion and absorption), exploiting a previously published glucose absorption formulation [40], which was demonstrated to adapt well to experimental data from individuals ranging from normal subjects to type-2 diabetic patients.

Sorensen model analysis
The Sorensen model [23] is one of the most complex physiological models describing glucose homeostasis. It has been vastly cited and used to simulate virtual patients with the aim of validating control algorithms in the framework of artificial pancreas development. It consists of three sub-models describing the glucose concentration time-course in brain, liver, heart and lungs, periphery (tissue and muscles), gut and kidney and includes the pancreatic release of insulin and glucagon. A scheme of the model is reported in Fig 1. The original model is represented by the first two block diagrams; the third diagram represents instead the added gastrointestinal tract, whose output constitutes one of the inputs into the G G Sorensen compartment. Appendix A1 reports the complete and corrected Sorensen model along with all the parameter values and the initial conditions of the equations. In presenting the model, Sorensen reported detailed reasoning and literature references for each equation and parameter value adopted, justifying each single choice with an accurate and broad discussion of the involved physiological aspects. A summary of the final model formulation was reported at page 213 of his dissertation in the subsection "Summary of Model Equations, Parameter Values, and Mathematical Nomenclature". Given its complexity, most of the researchers who used the model in their activity referred this summary section without double-checking the equations reported in the other sections of the original work. Some other Authors then referred to previous publications based on the Sorensen model, thereby inheriting the same imprecisions from their sources. Table 1 reports a list of errors we found in the

PLOS ONE
Glucose / Insulin Modelling original dissertation and which were incorporated in later work using the model [26,31,35,38,39]. Model simulations highlight the fact that the identified errors produce a clearly different model behavior: kidney glucose excretion appears to be slower (A), adopted initial conditions do not appear to be at equilibrium (C and E), insulin secretion appears to be incorrect (D). Error C was present for example in [26], [35], [38], [39] and in [31] where reference [35] was used as source. Error A was inherited by [26] and [38]. We hypothesize that the problem highlighted in this last publication, where equilibrium points appear to lie in a "non-feasible region", could in fact be due to the introduction of the above errors in the Sorensen model implementation used.

Sorensen re-implementation
Implementation and test of the original version. The Sorensen model was implemented by following the internal standard procedures of the CNR-IASI BioMatLab. Implementation of the model equations was conducted following the MoSpec (model specification) approach, an internally developed automatic system, which takes as input a spreadsheet containing all the model specifications (description of the equations, parameter names, initial conditions, parameter values, latex and matlab syntax) from which computational routines in Matlab [41], R [42] and C++ [43] are automatically generated, along with a LaTeX document containing the model equations as well as the program code produced for the implementation. Having automated code development and automatic comparison of the original mathematical formulation with the actual coding implementation allows fast and thorough code verification. In order to compare predictions obtained by Sorensen in his original work and predictions obtained by our system, a number of test cases presented by Sorensen were simulated: • a standard 0.5 g/kg Intra Venous Glucose Tolerance Test (IVGTT); • a variable-dose IVGTT comparison (0.05, 0.2, 0.5 and 0.75 g/kg); • a 0.04 U/kg Intravenous Insulin Tolerance Test IVITT; • a continuous intravenous insulin infusions (0.25, 0.4 mU/kg).
The next section reports the results obtained. The updated original version: Modelling the oral glucose assumption. The original work of Sorensen presents a series of simulations showing the good adaptation of the model to experimental data. Some of the simulations were those described above and show the time courses of glucose and insulin concentrations when glucose and/or insulin were administered intravenously as a bolus or as a constant infusion.
Sorensen also tested the performance of the model when glucose was administered orally: in this case he was forced to derive empirically both the rate of insulin secretion and the rate of glucose appearance. Two major limitations are indeed present in the Sorensen model: both the incretin effect and the rate of gastric uptake are not explicitly modeled. When glucose is ingested orally, glucose-induced insulin secretion by pancreatic β-cells is potentiated with respect to the case in which glucose is administered intravenously. This is due to the incretin effect of at least two hormones, glucose-dependent insulinotropic polypeptide (GIP) and glucagon-like peptide-1 (GLP-1), both released in the gut from enteroendocrine cells when glucose goes through the intestinal tract. This potentiation effect is not explicitly modelled by Sorensen and indeed the set of parameter values provided by Sorensen does not produce an insulin secretion rate compatible with the insulin concentrations that are normally observed during an oral glucose load. In order to obviate this discrepancy, Sorensen bypassed completely the pancreatic insulin release sub-model and iteratively adjusted the rate of assumed pancreatic insulin secretion into the circulating insulin sub-model, until the predicted serum insulin concentrations matched the experimental data. The same procedure was adopted for the rate of gut oral glucose absorption (r oga , mg/min), for which an empirical function was used as input in the glucose gut compartment and then adjusted so as to match the observed glucose concentrations with model predictions.
In the present work, an improved version of the Sorensen model is implemented, incorporating into the original 1978 version a previously published model of the Oral Glucose Tolerance Test, the SIMO model [40]. When choosing which representation of gastrointestinal glucose absorption to use in order to integrate the original Sorensen model, we reviewed several options from the literature. Appendix A2 describes the other gastrointestinal glucose absorption models that we considered. We finally choose the SIMO model due to its relative simplicity, to the fact that it represents mechanistically anatomical compartments and transfers, and to the fact that it fits available observations well. The SIMO model is composed of four compartments corresponding indicatively to Stomach, Jejunum and Ileum, plus a delay compartment between Jejunum and Ileum. Although the original SIMO model explicitly incorporates the incretin mechanism, making insulin secretion depend on the glucose content in the Jejunum and Ileum, the part of the model used here is only that related to the glucose absorption route. In the present formalization, stomach emptying and transfer from stomach to proximal and distal intestinal tract are modelled by linear dynamics. This approach was sufficient in the original SIMO model to obtain a very good adaptation of the model to data from subjects with different degree of glycemic impairment (from normal glucose tolerance to impaired fasting glucose and impaired glucose tolerance up to Type-2 Diabetes Mellitus). The equations inherited from the SIMO model are reported below: where D is the orally administered quantity of glucose whereas S, J and L represent the quantity of glucose in the stomach, in the jejunum and in the ileum compartments respectively. R is a delay compartment, necessary to approximate the absorption profile of glucose over time: see [40] and [44] for details. The function r oga represents the rate of gut oral glucose absorption, which acts as input into the original Sorensen model and which substitutes the empirical test function described above. In other words, the function r oga constitutes the link between the SIMO model and the Sorensen model (see Fig 1) and is substituted into the Sorensen model equation The incretin effect was not explicitely modelled by Sorensen and the original insulin model and pancreatic insulin release remain unchanged in the present formulation. Most of the parameters related to insulin metabolism are allow to vary freely, in order to describe a pancreatic insulin release able to produce a good fit to observed insulin concentrations. The next section reports the procedure carried out to test the ability of this improved version of the Sorensen model to reproduce an Oral Glucose Tolerance Test.

Model fitting and sensitivity analysis
The Sorensen insulin pancreatic release sub-model equations are reported below: The free parameters to be estimated are those reported in Eqs [1-5; 7-14], except for parameter Q 0 which was kept fixed at its original value as identified by Sorensen. For a detailed description of the Sorensen model variables and parameters refer to [23].
The parameter identification process was conducted according to four steps: The r oga time course, as empirically derived by Sorensen, was used as data input: data were graphically extracted (27 data points) and used to estimate the SIMO model parameters, so as to predict a curve that was as close as possible to the original Sorensen r oga function.

STEP 2
The SIMO model parameters were then kept fixed and one hundred optimization procedures were performed to estimate the free parameters of the Sorensen pancreatic insulin release sub-model, using as data glucose and insulin concentrations as well as points extracted from the r PIR time course (25 data points). In each one of the 100 optimization procedures the starting value used for each of the free parameter to be estimated was a random realization from a normal distribution centered on the parameter value reported by Sorensen with a standard deviation equal to 15% of the same parameter value.

STEP 3
A final optimization procedure was performed by letting all parameters from the above two steps vary freely, using all available data points from the four "observed" state variables (plasma glucose and insulin concentrations, r oga and r PIR ). This last optimization procedure was carried out using as starting value the optimum parameter vector from STEP 1 and the "best" optimum (the one producing the smallest loss function) among the 100 optimization procedures from STEP 2. An approximation to the Variance-Covariance matrix of the model parameter vector has been computed as [J T S −1 J] −1 , where J represents the Jacobian in the obtained optimum and S is Variance-Covariance matrix of the observation error vector, assumed to be diagonal (uncorrelated errors) with variances proportional to the square of the mean response: where σ, the scale parameter, is the coefficient of variation.

STEP 4
Two hundred optimization procedures were then performed to test model "a-posteriori" identifiability: each of the free parameters of the Sorensen insulin pancreatic release sub-model was initialized perturbing its optimum value from STEP 3: for each parameter to be estimated the new starting value in each of the 200 optimization procedures was sampled from a normal distribution centered on the parameter optimum value, with a standard deviation equal to its 15%.
Each optimization process was performed by minimizing the sum of the weighted squared residuals (weighted least-squares estimation, with weights the inverse of the squared expectations) over 5 free parameters in STEP 1 ( Table 2 with parameter f kept fixed at 1), 11 free parameters in STEP 2 (Table 3 with parameter Q 0 fixed at its original value) and 16 parameters in STEPs 3 and 4.
While it is true that the model incorporates a high number of free parameters to be simultaneously estimated and that this can lead to overfitting problems, STEPs 3 and 4 were executed in order to have reasonable starting values and reasonable allowable values ranges for each parameter (values for parameters can be assessed for "reasonableness" if each parameter has a direct physiological meaning, which is a good motivation for the use of mechanistic models). The results from STEP 4 give indications about both the variability of the parameter estimates and the correlation between each couple of parameter estimates. Table 2. STEP 1 results. Simo model parameters before (as reported in [40]) and after the optimization process for r oga data fitting.

Parameters [Units]
Before After

Bayesian a-posteriori identifiability analysis
For parameters with high coefficients of variation, obtained from the estimated asymptotic variance-covariance matrix of the model parameter vector (STEP 3), a bayesian approach was also implemented in order to evaluate their a-posteriori identifiability with an alternative approach to STEP 4, setting all the others parameters at their optimum value. A Monte Carlo Markov Chain (MCMC) approach was followed, using a Metropolis-Hastings algorithm [45] nested into the Gibbs sampling [46][47][48][49], making use of samples from the full conditional distributions. The Bayesian model specification starts with the assignment of a distributional form to the observation error vector, along with the specification of a variance-covariance structure. We assume normally distributed errors, with zero mean and the same variance-covariance matrix as the one adopted in the WLS approach: where y j is the observation at time t j , the errors � j , j = 1, . . ., n are assumed to be uncorrelated and where, in analogy with Eq 15, An a-prior inverse gamma distribution was assumed for the parameter σ 2 (σ 2 � GI( n 2 , t 2 )). The parameter θ was log-transformed and a prior uniform distribution was assigned to the parameter vector log(θ): p(log(θ)) � U D , where D is the region delimited by the intervals centred on the logarithm of the parameter optimum values and with limits equal to ±100% of the log-transformed values. The parameter vector, object of the bayesian inference, is therefore the parameter ξ = (η, σ 2 ) = (log(θ), σ 2 )).
The full conditional distribution for the scale parameter σ 2 is still an inverse gamma: with y being the observation vector, n the number of total observations and F the diagonal matrix with elements Due to non-linearity in the regression model, the full conditional distribution of η = log(θ) cannot be calculated explicitly. However it can be written up to a proportionality constant: The Metropolis-Hastings algorithm was used to obtain samples from the full conditional distribution of parameter the vector η = log(θ). At each iteration of the MCMC, a multinormal distribution for η was used as proposal distribution. Mean and variance-covariance matrix of the distribution were derived from the following considerations: using the first order Taylor expansion of the logarithm function, we can approximate log(x) in a small neighborhood of x, With the appropriate substitutions, the proposal multinormal distribution was assumed to be centred on the logarithm of the optimum values from the WLS procedure of STEP 3, while the variance-covariance matrix of η was derived from the coefficients of variation computed from the estimated asymptotic variance covariance matrix σ 2 [J T S −1 J] −1 . While the variance-covariance matrix remained unchanged during the MCMC algorithm, at each iteration of the algorithm the proposal distribution was centred on the values obtained from the preceding iteration. The MCMC algorithm was implemented by building nine independent chains of 10,000 samples, each one beginning from a different parameter starting point. The nine points were chosen as the four vertices of the region D plus the mid points of the edges and the central point of D. For each chain the first 1000 realizations were considered as "burn-in" to let the Markov chains to reach equilibrium and get sufficiently close to the stationary distribution. The characteristics of the empirical distribution of parameter θ can be derived from the characteristics of the empirical distribution of η. For each couple of parameters, a Credibility Region (CR) can be built from their bivariate empirical distribution: from a 20 × 20 bi-dimensional frequency histogram (400 cells in total with relative frequencies f i,j ), the frequency level ℓ such that the sum of the frequencies for which f i,j > ℓ is equal to 95%, is found iteratively: Deviations from the original time courses (Fig 73 at page 273 of [23]) are evident only for the insulin concentration profile, which is estimated to be about 20 mU/L lower for the 0.5g/ kg IVGTT experiment.  [23].
Here the glucose state variable is multiplied by 0.925 to obtain venous plasma glucose concentration, according to the prescriptions of the Author. The curves obtained again closely resemble, both quantitatively and qualitatively, the results by Sorensen.
Results related to the Continuous Intravenous Insulin Infusions experiment are reported in the bottom panel of Fig 3. In this last experiment, 0.25 mU/kg of insulin is infused intravenously for 150 minutes. Also in this case the state variable G PV (Glucose Periphery Vascular blood water space) is multiplied by 0.925 to obtained venous plasma glucose concentrations. From a comparison with Fig 75 at page 279 of [23], a very good superimposition of the obtained curves with the original ones is apparent. Table 2 reports the results derived from the first step of the identification process. Values in the second column are the average values reported in [40], which represent the starting point The top panel on the left reports the "observed" data and the predicted rate of appearance, showing the good adaptation of the SIMO model to the r oga values. Top panel on the right reports the "observed" (asterisks) and predicted values (continuous line) of pancreatic insulin release. The bottom panels report observed and predicted values of peripheral glucose (on the left) and insulin concentrations (on the right). The obtained predictions are very concentrated around their median curve and show a good adaptation of the model to data points. Table 3 reports the insulin release sub-model parameters of the new version of the Sorensen model (including the gastrointestinal tract) before (original values as reported in [23]) and after the optimization procedure from STEP 3. Estimates of many of the insulin secretion sub-model parameters differ from their original values, with a minimal difference of about 18.44% (β pir4 )  Table 3) and with the parameter values from STEP 3 (dashed line and continuous line respectively). Table 4 summarizes the 200 optimization procedures performed to check model identifiability: it reports the mean, standard deviation and coefficient of variations (CVs) of the parameter estimates obtained in the 200 procedures. CVs range from a minimum of 3.9% to a maximum of 22.2%.

Discussion
Much previous work addressed the need of developing glucose/insulin models, with different levels of complexity, for a variety of reasons, such as the study of insulin sensitivity or for controlled automatic insulin delivery (artificial pancreas) [7][8][9][10][11][12]. These models were to be identified on each single patient and had therefore to incorporate the relevant physiological mechanisms in a simplified fashion. More extended models, including a variety of interactions, are in principle more representative of the physiology, and have potentially a greater predictive ability. Such extended models, however, include a large number of parameters: most of them must be fixed to values taken from the literature in order for the remaining ones to be estimated from patient data.
As an example, the development of robust control algorithms for automatic glucose control, capable of maintaining glycemia within a normal range under different perturbations while minimizing the risk of dangerous hypoglicemic episodes, may make good use of adequate, realistic models. Motivated, among other things, by the need to provide plausible virtual diabetic patients for the development of efficient automatic glucose control laws, several complex, "extended" mathematical models of the glucose-insulin system have appeared in the last few years ( [17,20,23,26]).
The Sorensen model ( [23]) is perhaps the most detailed of these, in terms of documented parameter values and detailed physiological mechanisms. With its 22 nonlinear differential

PLOS ONE
equations and 135 parameters it has been vastly used in glucose control research. This model represents the time course of glucose concentrations in brain, liver, heart & lungs, periphery (tissue and muscles), gut and kidney, and includes submodels for the pancreatic release of insulin and glucagon. A more recent model, also belonging to the category of extended models, is the UVa/ Padova Type 1 Diabetes (T1D) Simulator ( [17][18][19][20]). Its last version incorporates some timevarying parameters, such as an insulin sensitivity function that changes in order to account for daily tissue glucose uptake variability or for the dawn phenomenon (observed elevated glucose concentrations in the very early hours of the morning). The UVa/Padova Simulator, in almost all of its several subsequent formulations, includes a glucose subsystem, a glucose rate of appearance subsystem, different routes for administered insulin, and a submodel for the glucagon kinetics and secretion.
The Hovorka model ( [21]), while substantially simpler than either the Sorensen or the UVa/Padova model, has also been used by several Authors for "closing the patient loop" in the development of automatic control laws [22,[50][51][52].
The rationale for using the Sorensen model, in particular for using an updated and corrected version of it, resides in the fact that it is a detailed, extended model for the simulation of

PLOS ONE
a virtual patient, that it is well documented (in terms of model equations and parameter values), and that is relatively easy to use it to replicate a wide range of experimental scenarios, as reported in Sorensen's original work.
Conversely, the UVa/Padova Simulator, which has been commercially distributed, is relatively difficult to implement independently because of the difficulty in finding what parameter

PLOS ONE
Glucose / Insulin Modelling values were used by the Authors in their simulations, and also because of the fact that some model equations are only partially described. The UVa/Padova is a simulator for Type I Diabetes Mellitus Patients, and for this reason it lacks any functional representation of the insulin secretion mechanism. In contrast, the Sorensen model allows the simulation of normal subjects, of type 2 Diabetes patients who still maintain some level of endogenous insulin secretion, and of Type I Diabetes patients for which endogenous insulin production is completely replaced by external inputs.
The Sorensen model however presents two major limitations: on one hand it lacks any representation of oral glucose administration and of the corresponding glucose rate of appearance in plasma; a second limitation is that insulin administration is foreseen only by means of intravenous boli or infusions. This latter is an important shortcoming: in clinical practice, subcutaneous delivery represents the normal route of insulin administration, by means of either subcutaneous boli or of wearable pumps. The subcutaneous insulin administration route could be introduced into the model as a further improvement.
In this work we begin with a thorough revision of the Sorensen model, obtaining a model implementation without the several imprecisions reported in the original work and inherited by the Authors who have used this model in their research activity. This corrected model implementation is made available to the scientific community in user-to-machine and machine-to-machine versions at the address: http://biomatlab.iasi.cnr.it/models/login.php (access as a Guest).
Matlab code is also downloadable from the same link. Moreover, with the aim of making the model more useful for researchers, we propose an improved version supplemented with a mathematical representation of gastrointestinal glucose absorption during oral administration. This updated version of the Sorensen model can also be found at above mentioned address. Following internal BioMatLab standards, the two model versions have been automatically implemented in Matlab, R and C++ starting from formal detailed instructions collected in a Model Specification (MoSpec) spreadsheet. The production documentation includes side-by-side computational code and mathematical description of the model equations, thus allowing a direct check of the correspondence of the computation with the intended mathematical formulation. In the final configuration, Matlab or R user interfaces (with good graphical capabilities and a potentially vast set of ready-made functions for post processing) exploit an underlying compiled fast C++ numerical integration engine. The compiled C++ engine also supports a Visualizer environment, where the user can quickly explore model behaviour corresponding to parameter changes.
The implementations in different languages and the Visualizer environment represent, collectively, a robust approach to verify the correctness of model implementation, which is the BioMatLab standard for model description, implementation and verification.
A series of simulations presented by Sorensen were replicated with our system in order to compare the two implementations. Simulations performed with the corrected version of the Sorensen model closely resemble the curves produced by the original model (see Figs 2 and 3), with very small deviations from what was reported in his work [23]. Small divergences are due to differences in some parameter values: some parameters used by Sorensen (such as basal glycemia or insulinemia) were said to derive from observed population average values but the actual numerical values adopted in the simulations were not reported.
In order to test the ability of the model, supplemented with the gastrointestinal tract submodel, to reproduce insulinemia and glycemia time courses from Oral Glucose Tolerance Tests (OGTTs), the improved version of the Sorensen model was tested on a set of OGTT data already presented by the Sorensen [23]. Fitting the improved model onto OGTT data shows a good qualitative behavior of the model solution with a good adaptation of the curves to experimental points (see Fig 4). We believe that two reasons are responsible for these good results: the first reason is that the Sorensen model, although built on a knowledge basis dating from the Eighties, is a well researched model, incorporating fundamental physiological mechanisms with reasonable mathematical representations. The second reason is that the added gastrointestinal model was taken from the SIMO publication ( [40]), where it was demonstrated that linear stomach emptying and linear glucose transfer to the intestine were sufficient to obtain a good representation of glycemia and insulinemia time courses data from subjects with different degree of glycemic impairment undergoing OGTT's [40].
Adapting the improved model to OGTT data required appropriate changes of some key parameter values of the Sorensen sub-model for insulin secretion. These changes were due to the fact that the Sorensen model does not include explicit modelling of the incretin effect [23]. Modifications of some parameter values however can compensate for this deficiency. It is well known ( [53][54][55][56][57]), in fact, that when administering glucose orally, the release of insulin by the pancreas is enhanced with respect to the same amount of glucose administered intravenously. This effect is due to the release of incretin hormones (GIP, GLP1) by the gut when it comes in contact with intraluminal glucose. The Sorensen insulin secretion sub-model does not contemplate an incretin effect, therefore, when using this model with its original parameter values, the model forecast of the insulin levels following oral glucose administration is lower than appropriate. On the other hand, since the glycemic increase after oral administration is slower than after intravenous administration, expected insulin release is also slower. In this case, a better approximation to the actual insulin levels can be obtained by allowing some of the parameter values to change; Table 3 reports the original parameter values and the values obtained after fitting the improved model onto OGTT data. For example, the parameter γ (Eq 10) which regulates glucose-driven transfer of insulin to the labile compartment is here estimated, in the case of an OGTT, to roughly half of the original value provided by Sorensen. This result is expected, and depends on the fact that when glucose is infused intravenously its absorption into the bloodstream occurs more rapidly, producing higher peak concentrations than when it is given orally. On the other hand, while γ is reduced, producing a slower increase of the releasable insulin Q, parameters β pir5 and β pir1 (Eqs 12 and 13 related to glucose driven early insulin release and to the potentiator factor respectively) are increased, producing a higher total release of the hormone.
The changes made to the parameter values determine changes in the dependent variables. For example, the total amount of labile insulin provision resulted to be increased (since for example β pir5 was increased). Changes in parameter β pir1 also determined large changes in glucose-enhanced early insulin release, X(G)): since the relationship is highly non-linear, an increment of about 30% of β pir1 results in a large increase of X(G). A higher release of insulin, in turn, requires a higher insulin inhibition in order to avoid hypoglycemias; this is represented in Sorensen's model by the dynamics of the inhibitor I (Eq 9), which is made to follow the insulin release X much faster by increasing the value of the parameter β.
The insulin secretion function (Eq 11) is governed by two terms: one acting at the early stage (M 2 (X − I)) and one acting at a later stage (M 1 Y). Optimization produced a smaller estimate for M 1 resulting in the function Y (the secretory effect of glucose on late insulin release, Eq 13) increasing more slowly for small values of G H (which occur at the beginning of the experiment), and more quickly when G H reaches higher concentrations (which, in OGTT's are maintained for a longer time in comparison with IVGTT's, requiring therefore an enhanced insulin production; see Eq 12). The reduction in the M 1 parameter may be necessary in order to balance the large increase in Y determined by the increase in β pir5 , particularly at high glycemias (at low glycemias the effects of changes in M 1 are negligible).
Moreover, a decrement of the parameter α in Eq 8, which represents the speed with which the potentiator factor P moves towards its target P 1 , contributes to maintain high level of insulin at later times.
Since the Sorensen model is very complex, even more so with the proposed addition of a gastrointestinal absorption submodel, it clearly presents identifiability problems. In fact, repeated fitting has shown that some parameter estimates are strongly correlated. In particular, parameters β pir1 and β pir4 present with a correlation coefficient of 0.99 (see Fig 6). This phenomenon is evident also by observing the results of the bayesian a-posteriori identifiability procedure, where the a-posteriori distribution of the two parameters that exhibit the largest coefficients of variation (β and M 2 ) is empirically derived. The 95% Credibility Region shows that in correspondence of small values of the parameter M 2 , the parameter β assumes indistinguishably small or large values. This is likely due to the fact that the action of the insulin inhibitor (which is expressed through the β parameter) exerts little influence for low insulin secretion (low values of parameter M 2 ), in particular during the initial stages of the experiments, where M 2 is most relevant, and where glucose concentrations have not yet reached high levels. Conversely, its contribution becomes more and more important as the secretion increases, so that an increase in M 2 makes an increase of β necessary. The descendent tract of the horseshoe corresponds to large values of β, when the Inhibitor equals the insulin release X(G) and the parameter M 2 has little influence (see Eq 11) and may take any value. We can see therefore how the bayesian approach helped to better identify the credible region for these two parameters.
It should finally be noticed that the Sorensen model makes use of hyperbolic tangent (tanh) functions, which are numerically effective, but which could be replaced by appropriate Hill functions, more immediately understandable by physiologists and clinicians.
In conclusion, in spite of its limits, the Sorensen model appears to be a flexible and well studied model, able to reproduce realistic physiological behavior from different experiments, particularly in the improved version discussed here. It is hoped that it will be found a useful instrument for the simulation of glucose/insulin system of virtual patients.

The Sorensen Model
Mass Balance-Glucose HEART AND LUNGS. GUT. LIVER. KIDNEY.
Metabolic Source and Sinks-Insulin Mass Balance-Glucagon Metabolic Source and Sinks-Glucagon

Initial Conditions-Glucose
Mass Balance. Metabolism.
ð101Þ P B ¼ P1 ð102Þ where: • q sto1 (t) and q sto2 (t) are the amounts of glucose in the stomach (solid phase and liquid phase, respectively); • δ(t) is the impulse function; • D is the amount of ingested glucose; • q gut is the glucose mass in the intestine; • k 21 is the rate of grinding; • k empt is the rate of gastric emptying; • k abs is the rate constant of glucose-intestinal absorption; • f is the fraction of glucose-intestinal absorption which appears in plasma; • k max is the maximum value for k empt ; • k min is the minimum value for k empt ; • α and β are the rates of decrease and increase, respectively, for the k empt function;

Elashoff model
The Elashoff model [59] assumed that the fraction of glucose in the duodenum compartment increases following a power exponential function. The equation system is as it follows: where: • q duo (t) is the glucose mass in the duodenum; • D is the amount of ingested glucose; • q gut is the glucose mass in the intestine; • k is the rate of emptying; • β is the shape factor.

Salinari et al. model
In the model reported in Salinari et al. [44], the transit of glucose through the intestine is represented by a mono-dimensional process where the glucose particles are transported from the proximal region to the distal region with a constant velocity. The system equation is the following: @q @t þ u @q where: • q is the glucose density; • γ is the rate of glucose absorption.
The gastric emptying is given by the boundary condition in z = 0: where: • ZðtÞ ¼ Dbk b t bÀ 1 e À ðktÞ b is the rate of glucose delivery in the duodenum; • D is the glucose dose; • u is the constant velocity with which the glucose particles are transported from the proximal region to the distal region.

Lehmann & Deutsch model
The Lehmann & Deutsch model [60] describes glucose absorption by the gut, assuming that the gastric emptying process is represented by a trapezoidal function. The absorption in the intestine follows a first order kinetics. The system equations are: where • q gut (t) is the amount of glucose in the gut; • k abs is the rate constant of glucose intestinal absorption; • f is the fraction of the glucose intestinal absorption which appears in plasma; • Ra is the rate of glucose absorption.