A modelling approach to hepatic glucose production estimation

Stable isotopes are currently used to measure glucose fluxes responsible for observed glucose concentrations, providing information on hepatic and peripheral insulin sensitivity. The determination of glucose turnover, along with fasting and postprandial glucose concentrations, is relevant for inferring insulin sensitivity levels. At equilibrium (e.g. during the fasting state) the rate of glucose entering the circulation equals its rate of disappearance from the circulation. If under these conditions tracer is infused at a constant rate and Specific Activity (SA) or Tracer to Tracee (TTR) ratio is computed, the Rate of Appearance (RA) equals the Rate of Disappearance (RD) and equals the ratio between infusion rate and TTR or SA. In the post-prandial situation or during perturbation studies, however, estimation of RA and RD becomes more complex because they are not necessarily equal and, furthermore, may vary over time due to gastric emptying, glucose absorption, appearance of ingested or infused glucose, variations of EGP and glucose disappearance. Up to now, the most commonly used approach to compute RA, RD and EGP has been the single-pool model by Steele. Several authors, however, report pitfalls in the use of this method, such as “paradoxical” increase in EGP immediately after meal ingestion and “negative” rates of EGP. Different attempts have been made to reduce the impact of these errors, but the same problems are still encountered. In the present work a completely different approach is proposed, where cold and labeled [6, 6-2H2] glucose observations are simultaneously fitted and where both RD and EGP are represented by simple but reasonable functions. As an example, this approach is applied to an intra-venous experiment, where cold glucose is infused at variable rates to reproduce a desired glycaemic time-course. The goal of the present work is to show that appropriate, if simple, modelling of the whole infusion procedure together with the underlying physiological system allows robust estimation of EGP with single-tracer administration, without the artefacts produced by the Steele method.


Introduction
Insulin resistance is the common denominator of many pathological conditions including type 2 diabetes (T2DM), obesity, hypertension, hyperlipidemia, atherosclerosis, non-alcoholic fatty liver disease (NAFLD) and polycystic ovarian syndrome. Drugs that improve insulin sensitivity, such as metformin or thiazolidinediones, improve glycemic control by acting in different ways. For instance, while metformin suppresses Endogenous Glucose Production (EGP), thiazolidinediones mainly increase peripheral glucose uptake [1]. Bariatric surgery (of different types) determines diabetes remission both in the short [2,3] and in the long term [4][5][6] through reversal of insulin resistance [7]. However, while Roux-en-Y Gastric Bypass mainly reduces hepatic insulin resistance, Bilio-Pancreatic Diversion increases whole-body insulin sensitivity [8,9], which is mostly determined by skeletal muscle mass [10]. In this context it is clear that the assessment of the differential mechanisms, whereby insulin sensitivity is improved, greatly benefits from the possibility of estimating Endogenous Glucose Production (EGP). Stable isotopes are currently used to measure glucose fluxes responsible for observed glucose concentrations, providing information on hepatic and peripheral insulin sensitivity: in particular, hepatic insulin sensitivity can be gauged by measuring EGP. During fasting, EGP provides the necessary glucose to maintain euglycemia: the liver, through glycogenolysis and gluconeogenesis, contributes to about 90% of EGP [11], while kidneys and gut are jointly responsible for the remaining 10%. In the postprandial state the liver helps maintain normal glucose levels by removing glucose from the circulatory stream and storing it into glycogen, while EGP is suppressed. For constant plasma insulin, glucagon and growth hormone concentrations, a doubling of plasma glucose levels was observed to inhibit EGP by 42% and stimulate peripheral glucose uptake by 69% in nondiabetic subjects; in T2DM patients EGP was not suppressed and glucose uptake was stimulated by only 49% [12]. Portal venous hyperglycemia and hyperinsulinemia stimulate hepatic glucose uptake [13] and promote net hepatic glycogen synthesis [14]. Experimental data have indicated that extrahepatic signals may also be important in influencing EGP. Lactate and alanine [15], but also glycerol and free fatty acids [16], are gluconeogenetic precursors: higher concentrations of these substrates determine an increase in EGP. Insulin inhibits gluconeogenesis while glucagon and catecholamines stimulate it [17]. In T2DM hyperglucagonemia contributes to EGP and hyperglycemia [18].
The determination of glucose turnover, along with fasting and postprandial glucose concentrations, is relevant for inferring insulin sensitivity levels. At equilibrium (e.g. during the fasting state) the rate of glucose entering the circulation equals its rate of disappearance from the circulation. If under these conditions tracer is infused at a constant rate and Specific Activity (SA) or Tracer to Tracee (TTR) ratio is computed, the Rate of Appearance (RA) equals the Rate of Disappearance (RD) and equals the ratio between infusion rate and TTR or SA. In the post-prandial situation or during perturbation studies, however, estimation of RA and RD becomes more complex because they are not necessary equal and they may also vary over time due to gastric emptying, glucose absorption, appearance of ingested or infused glucose, variations of EGP and glucose disappearance.
Up to the present day, the most commonly used approach to compute RA, RD and EGP has been the single-pool model by Steele [19][20][21], which has been used by many investigators [22][23][24][25]. However, several authors report pitfalls in the use of this method, such as the "paradoxical" increase in EGP immediately after meal ingestion and, more troubling, "negative" rates of EGP [26,27].
One may attempt to reduce the impact of these errors by using experimental designs where SA or TTR are kept as constant as possible and where sampling is very frequent [28], but this would limit the range of practicable designs and red would increases experimental costs. Other authors have proposed different approaches: Radziuk et al. advocated the use of a two-compartment model [29], with the introduction of an accessible and a slowly exchanging nonaccessible compartment, but the same problems were encountered as with the dual-tracer approach. Cobelli et al. [26] proposed the use of three tracers for the estimation of RA meal , EGP and RD after the ingestion of a meal: this method is however very complicated in that it requires, in addition to the use of a third tracer, also that pilot studies are conducted and that the intra-venous tracer infusion rates are varied over time so as to keep both oral and endogenous TTR's constant.
In the present work a completely different, novel model-based approach is proposed, where cold and labeled [6, 6-2H2] glucose observations are simultaneously fitted and where both RD as well as EGP are represented by simple but reasonable functions. As an example, this approach is applied to an intra-venous experiment, where cold glucose is infused at variable rates to reproduce a desired glycemic time-course. The goal of the present work is to show that appropriate, if simple, modelling of the whole infusion procedure together with the underlying physiological system allows for a robust estimation of EGP with a single-tracer administration (as in the Steele approach), without the additional complications of double or triple tracer administration and without the artifacts produced by the Steele method itself.

Study sample
Two NGT, two IGT and two T2DM subjects (average BMI of 53.5 kg/m 2 ) participated in an isoglycemic experiment according to a protocol (Trial registration: ClinicalTrials.gov NCT03223129) approved by the Ethical Commitee of the Catholic University School of Medicine Policlinico Gemelli, Rome, Italy. All participants provided written informed consent. Anthropometric characteristics of the six subjects are reported in Table 1. Details of the experiment and procedures are reported in [30]. Briefly, the experiment consisted in administering patients an oral and an intra-venous glucose test. All subjects were studied in the post-absorptive state after a 12-h overnight fast. To collect arterialized venous blood, a retrograde catheter was inserted in a dorsal hand vein, with the hand kept in a warming blanket. A forearm vein of the controlateral arm was catheterized for the infusions. On day 1, At 08,00 h, [6, 6-2H2]glucose was infused as a priming dose of 22 μmol/kg and then as continuous infusion throughout the whole procedure at 0.22 μmol/kg/min. After 2.5 h of isotope infusion (basal period), oral glucose loads were given and consumed over 5 minutes at intervals of two hours, three times. The three oral loads consisted respectively of 25 g, 75 g and 100 g of glucose in aqueous solution. Glucose concentrations were evaluated before the start of the first oral glucose load and every ten minutes thereafter. On day 2, the procedure started again, similarly to the first day, with [6, 6-2H2] glucose infusion, first with a priming dose of 22 μmol/kg and then as a continuous infusion throughout the whole procedure at 0.22 μmol/kg/min. After 2.5 h of [6, 6-2H2] glucose continuous infusion (basal period), a 20%wt/vol adjustable glucose solution was infused at varying rates to match the plasma glucose concentrations obtained during the oral glucose tolerance test performed on the previous day. The 20% dextrose that was infused during the intra-venous glucose infusion test was enriched with 2, 5% of [6, 6-2H2]glucose to minimize changes in glucose isotopic enrichment. During this experimental IV procedure plasma glucose was measured before the start of the adjustable glucose infusion and every five minutes thereafter in order to change the glucose infusion rate so as to obtain the "isoglycemic" pattern. Plasma glucose concentrations were determined by the glucose oxidase method.

Modelling
The objective of the present work was that of developing a model-based approach to estimate glucose liver production during the intravenous experiment described above or similar ones. According to the experimental procedure described in the previous paragraph, the data analysed in the present work is that collected on day 2 when subjects underwent the intra-venous glucose infusion reproducing the glycaemic Oral Glucose Tolerance Test (OGTT) patterns derived from the administration of the three oral glucose loads on day 1 (see Fig 1). The proposed modelling approach, therefore, is applied to the only intra-venous infusion procedure which, from a modelling point of view, is the simplest part of the whole experiment. References

PLOS ONE
to OGTTs are made, therefore, to highlight the three different trends (mirroring the three OGTTs) over time. We note that the proposed modelling approach is independent from the pattern of the infusion, and that, independently from the OGTTs performed on the previous day, it aims at modelling an intra-venous glucose experiment, by implementing the intravenous inputs, and the consequent physiological response. We start by considering that in humans water represents about 65% of total body weight (BW) and approximately 5% of BW is water in plasma, 20% in interstitial fluid and 40% in intra-cellular fluid. For a 70 kg man we can therefore state that about 3.5 ℓ plus 14 ℓ represent the distribution volume for glucose, about 0.25 ℓ/kg overall. In obese subjects the distribution volume per kilogram body weight is generally smaller.
An average consumption of glucose by the brain of about 104 g/day is considered [31,32]. This means about 0.40 mmol/min (104/1440 × 1000/180 mmol/min). Upon administration of tracer, the total glucose uptake (sum of 'cold' and 'hot' glucose uptake, where with 'hot' we indicate the [6, 6-2H2]glucose, and with 'cold' we indicate the tracee) by the brain can be divided into two components. The hot component, k brh,fun /V gp in equations (4) is the proportion TTR/(1+TTR), where TTR represents the tracer/tracee ratio: k brh;fun ðtÞ=V gp ¼ ð104=1440 � 1000=180Þ � TTRðtÞ=ð1 þ TTRðtÞÞ which varies over time as function of TTR(t), while the total uptake in Eq (1), represented by the term k br /V gp is constant with k br equal to 0.40 mmol/min (being V gp the plasma glucose distribution volume expressed in ℓ/kg).
Two different models were evaluated: the simplest model includes only one physical glucose compartment, while a more detailed model includes both plasma and interstitial compartments. When evaluating the last model two volumes were thus separately considered: plasma volume, whose typical value is approximately 0.05 ℓ/kgBW, and interstitial volume, whose typical value is about 0.20 ℓ/kgBW.
All the model state variables and model parameters are described in Tables 2 and 3, respectively, along with their corresponding units of measurement.
At time -150 min, corresponding to the time of the priming infusion, the system is supposed to be at steady state and all the starting conditions of the model equations are referred to  'Total' glucose concentration in plasma. The variation of 'total' glucose over time ('hot' plus 'cold' glucose and indicated with G P (t) in the equation below) is determined by the external glucose infusion (R inf ), by the insulin-dependent uptake of glucose by tissues (R upt , function of G P (t) and of the plasma insulin concentration I(t)), by the insulin-independent brain uptake (k br ) and by the term related to the exchanges between plasma and interstitium when two glucose compartments are considered (with G I (t) the total glucose concentration in the interstitium compartment). Because we consider experiments where very high glucose concentrations are not reached, we have set renal glucose extraction (represented in the model with the term R kid ) to zero. While in principle this choice could be simplistic, plasma glucose concentration in fact never exceeded the "threshold" of 11 mM [33] during the analyzed experiments. Moreover, because in a more general representation of the glucose dynamic the rate of appearance from the gastrointestinal tract could also be considered, this has been introduced in the formalization and indicated with R a,gut , set to zero here because no oral administration occurs.
In the presence of additional glucose input from the gastrointestinal tract, or as part of modeling a simple OGTT experiment, R a,gut can be modeled as f�G gut (t)/V g p. Here f is the fraction of the amount of glucose administered that is effectively absorbed through the intestine, and G gut (t) is the glucose that goes through the intestinal tract, possibly approximated as a series of consecutive compartments representing the length of the entire intestine, starting from the proximal down to the most distal part. Endogenous glucose production is represented by the term R a,liv in Eq (1). The expression for 'total' glucose variation is therefore the following: where R upt [mM] is a piecewise function: with R upt ðG P ðÀ 150Þ; IðÀ 150ÞÞ ¼ R uptb k xGI;1 � G P ðÀ 150Þ � IðÀ 150Þ: The function in Eq (2) depends upon three different constant insulin-dependent glucose elimination rates k xGI,j (j = 1,2,3), one for each of the three experimental sub-periods. We hypothesize, therefore, that the insulin sensitivity of a subject may vary over the course of the whole procedure: glucose administration, causing insulin release, would thus have different effects when given repeatedly. While some hypoglycemia-preventing mechanism, limiting over time the effect of circulating insulin, may be associated with oral glucose administration, the actual mechanisms involved in the variation of insulin sensitivity are not the object of the present work, but deserve to be investigated.
The terms −k ig G P (t) and k gi V gi V gp G I ðtÞ in Eq (1) represent the exchanges between the plasma and the interstitium compartment, with parameter k ig indicating the rate of transfer from the plasma to the interstitium and parameter k gi the constant transfer rate from the interstitium to plasma.
Finally, in the one-compartment glucose model, since there is no interstitial compartment, the last two terms in Eq (1) must be eliminated.
'Total' glucose concentration in the interstitium. Different mathematical models for the description of the equilibrium between plasma glucose versus interstitial glucose have been proposed in literature. These models are based on the assumption of free diffusion of glucose molecules between plasma and interstitium, and of glucose uptake from the interstitial fluid by tissue cells [34][35][36][37]. Rebrin and Steil proposed a "two-compartment" formulation, where interstitium and plasma are two independent compartments separated by a diffusional barrier which glucose is free to traverse, based on its concentration gradient [36,38,39]. In this formulation glucose is cleared from the interstitial fluid by tissue cells, at a rate dependent on the interstitial glucose concentration itself. We adopted essentially the same formulation, assuming no external loss from the interstitium. Moreover we hypothesized that the two compartments are at equilibrium at steady state so that their concentrations before the beginning of the experimental procedure are the same. Transfer of glucose from plasma to interstitium and vice versa is described by the transfer rate constants k ig and k gi .
From equilibrium conditions on Eqs (1) and (3) it follows that two parameters can be determined: where R ab,liv is the basal endogenous glucose production. The last equation is also valid for the one-compartment glucose model. Plasma concentration of infused 'hot' [6, 6− 2 H 2 ] glucose. The isotopic hypothesis is that 'hot' [6, 6− 2 H 2 ] glucose has the same metabolic destiny as 'cold' glucose, so that, apart from the terms related to the Endogenous Glucose Production and forcing input functions, the variations of 'hot' glucose can be described, in analogy with Eq (1), by the following equation: where As described before, a priming dose of 22 μmol/kg was infused at the start of the experimental procedure (t = −150 min) and it is represented in the model as the initial condition in Eq (4); the term R inf,H represents instead the continuous infusion of 0.22 μmol/kg −1 administered throughout the whole procedure, as described in the subsection "Study Sample".
Interstitial concentration of infused [6, 6− 2 H 2 ] glucose. In the same way as for 'total' glucose, a linear transfer between plasma and interstitial compartments is considered for 'hot' glucose, with the difference that while 'total' glucose is assumed at equilibrium at time t 0 , 'hot' glucose starts at t 0 from a concentration of zero, introducing a delay with which 'hot' interstitial glucose concentrations follow 'hot' plasma glucose concentrations: Tracer ([6, 6− 2 H 2 ] glucose) to tracee ratio. An external substance, the tracer, is used to follow a natural substrate, the tracee, when such a substance exhibits the identical metabolic fate as the tracee. Tracers used in human research are usually identical to the tracee, except that one or more atoms differ from the natural form. The assumption that tracers, isotopes, are treated in identical fashion as their counterparts is at the basis of tracer methodology. The isotopes fall into two classes: radioactive and stable. When dealing with radioactive isotopes the unit of enrichment is the specific activity, SA: Since scintillation counting is able to detect small quantities of radioisotope, the infused tracer is essentially "massless" and the amount of tracer is represented by disintegrations per minute (dpm) [40]. When dealing with stable isotopes, large amounts of labeled compound are necessary. The analogue of SA when stable isotopes are considered is the Tracer to Tracee Ratio (TTR), that is the ratio of the concentrations of tracer to tracee. From time 0 onwards, the Tracer to Tracee ratio (TTR) is measured every 20 minutes and in the present modelling is given by the ratio of 'hot' glucose over 'cold' glucose: Y represents the model prediction of the observed TTR's. Because 'cold' glucose concentrations were measured every 10 minutes, missing data were estimated by interpolation. From observed TTR's and observed glucose concentrations, observed 'hot' glucose concentrations were derived and used in the parameter estimation procedure. Estimation was performed by minimizing a loss function incorporating deviations between observed 'total' glucose concentrations, 'hot' glucose concentrations and their respective model predictions, see below.
Infusion of 'cold' glucose. Infusion of 'cold' glucose appears in both Eqs (1) and (4) and is indicated with R inf . It is not properly only 'cold' glucose since the 20% dextrose infused was enriched with 0,25% of [6, 6-2H2] glucose. This forcing input function is derived from the data, transformed to match the unit of measurements adopted by the experimenters and those adopted in the present work (i.e. from ml/h to mmol/kg/min). Considering that the infused solution was at 20%wt/vol, measurements were multiplied by the factor M inf = 1000×0.2/ (180×Weight × 60), where Weight is the subject's Body Weight in kg.
Rate of liver glucose production. Two formulations of the Rate of Appearance from the liver (that is Endogenous Glucose Production, EGP) have been tested. Liver glucose production from precursors is mediated by the inhibitory effects of insulin and of glucose itself. Hyperinsulinemia is moreover accompanied by suppression of glycogenolysis. Even if some studies suggest that insulin acts on EGP suppression through an extrahepatic route, such as via the suppression of the liberation of FFA's from adipose tissue [41,42], the adopted formulations synthetically model the effect of insulin and glucose directly on EGP.
A first mathematical formalization may adopt multiplication effects: EGP decreases from a maximal production rate, indicated with k IG,max ([mmol/kg/min]), towards zero for increasing insulin and glucose concentrations according to the following equation: which, from the steady-state condition derived from Eq (1), leads to the following expression for the parameter k IG,max : Adopting the above formulation requires the estimation of four free parameters related to EGP modeling. With the aim of reducing the number of free parameters to be estimated, so as to reduce identifiability problems due to overparametrization, a simpler formulation, requiring only two parameters, has been evaluated and then adopted as the definitive formulation: which, from the steady state condition of Eq (1) implies: Given that the parameter k IG,max can be thus determined, the single free parameter to be estimated in Eq (8) is λ GI .
Infusion of 'hot' glucose. After a priming dose of [6, 6-2H2]glucose, administered at time t=-150, a continuous infusion of 'hot' glucose is given during the whole procedure. The forcing function in this modelling formulation has been indicated with R inf,H ([mmol/kg/min]) and appears as the first term in Eq (4), constant throughout the whole procedure. The 0.25% of the forcing function R inf contributes to the 'hot' glucose infusion. The priming dose has been introduced instead in the model as the initial condition of the same equation.
Insulin plasma concentration. Insulin concentrations were measured at times 0, 20, 60, 100, 140, 180, 220, 260, 300, 340 minutes. Modelling of insulin secretion would in all likelihood prove beneficial from the point of view of overall model robustness, but a detailed discussion of this aspect would cloud the present issue. For the purpose of this work we have therefore performed a linear interpolation of the observed measurements, smoothing the obtained vector with a Gaussian-weighted moving average filter. The insulin forcing function was denoted with I(t) ([pM]).

Estimation
Model parameters, their descriptions, units of measurements and corresponding plausible values are reported in Table 3. The last column of the table shows if the parameter is "free" (and hence to be estimated), "determined" from the steady state conditions or "computed" on the basis of initial conditions. When the parameter is "determined", the corresponding equation in the manuscript is also reported. When it is "computed", the computation is shown in the table . Fig 1 shows a schematic representation of both the experimental procedure and its modelling approach for the one-compartment glucose model. The model was adapted to data by minimizing the weighted sum of squared deviations of the observed 'total' glucose concentrations and 'hot' glucose concentrations data from their respective predictions, with weights the inverse of the squared predictions. Estimation was also performed introducing serial correlation within each individual. The hypothesised correlation matrix is: where observations y i1 and y i2 on the i−th subject are taken at times t i,1 and t i,2 , respectively. Such a matrix accommodates situations of not equally spaced observations and of a correlation that decreases with increasing lags. If the model for correlation is given by (10) then the covariance for the observation vector y i is specified by the following expression: with θ the parameter vector of the model free parameters. Estimation of parameters σ, α and θ was carried out by a recursive procedure which implements subsequent steps: 1. derive an initial estimatesŷ of θ by means of an Ordinary Least Squares or Weighted Least Squares approach; 2. useŷ to obtain an estimate for α and σ with the pseudolikelihood method which minimizes the following function: 3. compute the covariance matrix Cov withâ andŝ from the previous step in (11) and use it to perform a Weighted Least Squares estimation with weights in Cov −1 4. repeat step 2) and step 3) until convergence.
In any case the final estimation of σ 2 is computed on the estimation of θ and α according: Under the hypothesis of normally distributed errors, the approximate covariance matrix Sŷ of y is given by: where J is the Jacobian, that is the (n×p) matrix with jth row equal to f ðx j ; yÞ 0 , with n and p the number of observations and the number of free parameters, respectively.
Free parameters to be estimated in the final adopted mathematical formulation were V gp , k xGI,j , for j = 1, 2, 3, G Pb and λ GI . When using the two-compartments glucose model, two more parameters had to be estimated: V gi and k ig . The use of the multiplication effect would required the estimation of three additional parameters: four new parameters would appear (K I , K G , γ I and γ G ), while one (λ GI ) would disappear. The remaining parameters were determined by steady state conditions or were set to values either known from the literature (this is the case for the glucose uptake by the brain, see the "Modelling" section above) or derived according to the experimental design (such as some initial conditions).

Results
The results related to the model including multiplicative effects have not been reported: although the fittings obtained were good, parameter estimates exhibited large variations among the studied subjects, and good model performance on the same subject could be obtained with very different values of the four parameters (K G , K I , γ G , γ I ), leading to the conclusion that this model was overparametrized. All results detailed below refer therefore to the model without multiplicative effects. Fig 2 shows the EGP computation from the Steele equation, both in the case TTR data are smoothed and in the case they are only interpolated (continuous green line and red dashed line, respectively, in Panel A). The figure reports four different trends of the EGP calculated by also varying the volume of distribution for glucose. According to [43], reasonable values for the volume could range indeed from 40 (plasma volume) to 230 ml/kg (extracellular fluid volume). The value often used is 145 ml/kg, which corresponds to a proportion p = 0.63 of the total volume (i.e., 230 ml/kg). The figure shows the EGP trend obtained with raw data after they have been only interpolated in the entire interval of observation (dashed red line) by using the standard volume of 145 ml; the dashed blue and black lines refer to computations performed by using smoothed interpolated TTR data, as well as smoothed Ra and EGP data, with distribution volumes equal to 145 ml and 40 ml, respectively; the continuous black line represents the predicted EGP estimated with the proposed modelling approach.
Figs 3-8 report the observations and the predictions of 'total' and 'hot' plasma glucose concentrations, as well as of TTR, along with the model predicted EGP (continuous lines)

PLOS ONE
superimposed to the estimated EGP as derived by the Steele equation (dashed line). The last two panels report smoothed insulin concentrations and the estimated rate of glucose disappearance (given by the sum of insulin-dependent glucose uptake and by the zero-order brain uptake) for the six studied subjects under the one-compartment glucose model. While it is

PLOS ONE
evident the good performance of the model in the studied subjects, it is especially important to consider the predictions of EGP. The mathematical formulation adopted produces predictions that can only assume positive values; furthermore, the paradoxical increase in EGP in the first minutes after the start of the experiment never occurs with the adopted model formulations.

PLOS ONE
Panels from A to E show two different prediction lines, black and blue, derived from WLS parameter estimation procedure and from the approach taking into account autocorrelated observation errors. When the two procedures produce very similar results, the two curves in each panel are superimposed and only the blue line is evident.
Of interest is the comparison between the performance of the two model formulations including or excluding the interstitial compartment. Figs 9 through 14 show the plasma glucose G P , the 'hot' glucose (H P ) and TTR prediction together with their respective observations, as well as the estimated EGP for all patients with model A (not including the interstitial compartment, right panels) and model B (including the interstitial compartment, left panels). Table 4 reports the estimated distribution volumes (only plasmatic, V gpA , for the one-compartment glucose model and both plasmatic and interstitial, V gpB and V giB , for the two-compartments glucose model) along with the transfer rate constants from the two compartments (plasma and interstitium) for model B as well as the loss function values, the Akaike (AIC) and BIC criteria obtained with the two formulations.
Both the AIC and BIC values privileged the simpler model (Table 4) and from the visual inspection (Figs 9-14) it is clear that in general little gain is obtained passing from the onecompartment glucose to the more complex two-compartment glucose model, except for Subject 2 for which the simpler model seems to perform better (Fig 10, panels A and C for the two-compartment glucose model and panels B and D for the one-compartment glucose model). Table 5 reports the estimates of the free parameters of the one-compartment glucose model under the WLS approach, whereas Table 6 reports the parameter estimates under the hypothesis of correlated errors. Fig 15 shows the scatter plot of weighted residuals towards observations. The average estimates and standard deviations derived from 14 are: V gp = 0.246±0.036; k xGI,1 = 3.67×10 −05 ±1.42×10 −05 ; k xGI,2 = 2.37×10 −05 ±9.42×10 −06 ; k xGI,3 = 4.09×10 −05 ±1.13×10 −05 ; λ GI = 0.0027±0.0013; G Pb = 6.10±3.24. Values of the three indices of insulin sensitivity highlight for each studied subject a decreasing trend from the first to the second phase of the experimental procedure and an increasing trend from the second to the third phase, suggesting an intraday variability of the insulin sensitivity. On the contrary, Figs 4-7 clearly show the occurrence of the inappropriate initial EGP peak when using the Steele equation, which also produces in all subjects the typical (erroneous) negative EGP predictions. The EGP time courses obtained with the modelling approach proposed in the present work are in accordance with the observed and predicted glucose concentrations: after an initial fall from the basal condition (due to the start of the adjustable glucose infusion), the model-predicted EGP remains mostly constant or presents slight oscillations for the whole experimental procedure, increasing over the intervals where glucose concentrations are restored to their basal values. It is interesting, in this respect, to observe that in Fig 7 a persistent increase in glucose concentration corresponds to nearly zero liver glucose production.

Table 4. Distribution volumes, loss, Akaike (AIC) and BIC criterium in the two model formulations: one-glucose compartment model (A) vs glucose two-compartments model (B).
For Model B also the transfer rates between plasma and interstitium are reported.

Discussion
Determination of glucose turnover, fasting and postprandial glucose concentrations is important for the assessment of possible impairment of glucose/insulin metabolism. Insulin resistance, both central and peripheral, plays an important role in the development of type 2 diabetes [44][45][46]. The liver is crucial for the maintenance of normal glucose homeostasis; it produces glucose during fasting and stores glucose as glycogen in the postprandial state: net hepatic glucose production derives from the sum of glucose fluxes determined by gluconeogenesis, glycogenolysis, glycogen synthesis, glycolysis and other pathways. In the postprandial state the liver helps to maintain normal glucose tolerance by absorbing glucose from plasma and storing it into glycogen, while hepatic glucose production (EGP) is suppressed. Stable isotopes are currently used to measure glucose fluxes and for providing important information on hepatic and peripheral insulin resistance. Tracer techniques are commonly employed to estimate the rate of Appearance (RA) and the Rate of Disappearance (RD). In the post-prandial situation or during perturbation studies, estimation of RA and RD is rather complex due to gastric emptying, glucose absorption, appearance of ingested or infused glucose, variations of EGP and of glucose uptake.
The single-pool model by Steele [19][20][21] has been for a long time the most commonly employed approach for estimating RA, RD and EGP and is still widely employed in clinical investigations. However, the method suffers from at least two problems: on one hand it produces "paradoxical" increases in EGP immediately after glucose intake or after a meal, when, on the contrary, physiology would predict a compensatory fall in EGP. The Steele approach also leads to the computation of "negative" rates of EGP: it is to be kept in mind, in this context, that the EGP we refer to here is not a net glucose production balance by the liver, but a true hepatic glucose output, and that therefore negative EGP values are by definition not

PLOS ONE
possible. Experimental designs with frequent sampling [28] have been employed to reduce these errors; Radziuk [29] suggested the use of a two-compartment model incorporating an accessible and a slowly exchanging nonaccessible compartment; other authors followed the dual-tracer approach, but the same problems were still encountered. Cobelli et al. [26] even proposed a more complicated method involving the use of three tracers for the estimation of the Rate of appearance (Ra meal ), EGP and R D after the ingestion of a meal together with variable IV tracer infusion rates so as to keep both all TTR's constant. All of these methods are more expensive in terms of resources and human labor, and they do not satisfactorily solve the original problems anyway.
With the object of overcoming such problems, in this work a completely model-based approach is followed, using a simple compartmental model for glucose (either including or not including an explicit interstitial compartment). The proposed model is actually only one out of many possible models that can be used for the analysis of stable isotope or radioactive data, with the aim of estimating endogenous glucose production. What is here advocated is not the use of this particular model (even though this model appears to work well, at least for the present experimental application), but rather a shift to an altogether different approach: the simultaneous estimation of both 'total' and 'hot' glucose concentrations produces a more accurate assessment of tissue glucose uptake and of EGP, by representing them directly with simple but reasonable functions. By considering jointly the time courses of 'hot' and 'total' glucose concentrations, the glucose rate of disappearance can be estimated more precisely from the dynamics of 'hot' glucose, since RA is the only unknown in the equation. The RA of 'cold' glucose is then derived directly from the isotopic assumption, that 'cold' and 'hot' glucose are subject to the same metabolic processes.
To exemplify this approach, we applied it to an intra-venous experiment, where 'cold' glucose is infused at variable rates to reproduce a desired glycemic time-course: in this experiment, we have to deal with an intrinsically non-steady-state situation, since the glycemic profile we are attempting to reproduce is itself not constant.
We tested two different mathematical formulations: one including only one compartment for glucose distribution in the body (the plasma compartment) and one including an additional compartment for glucose distribution in the interstitial space. The one-compartment glucose model has only six free parameters to be identified from glucose concentration and TTR data, considering insulin as a forcing function. The introduction of the interstitial compartment with which plasma glucose reaches equilibrium through free diffusion (determined by the concentration gradient as formulated also in Rebrin and Steil [38]) introduces two more parameters to be estimated.
Judging on the basis of the present dataset and the currently estimated parameter values it is not clear whether the additional compartment is necessary. The results highlight a significant increase in fitting performance with the two-compartment glucose model only for patient 1 (see loss values in Table 4), with an equivalent EGP time course (see Fig 9). In two cases, in the face of an essentially equal performance of the two formulations (see for example Figs 11 and 12 for Subject 3 and 4 respectively), the estimated parameters suggest that the additional compartment may produce not plausible parameter values: the optimization procedure leads to smaller values for the interstitial distribution volume V gi than expected, and produces smaller values for V gi than for the plasmatic volume V gp . Moreover, for Subject 2 the simpler model produces a better data interpretation (as shown in Fig 4 and in Table 4). Given these considerations, along with the computed AIC and BIC values, the present observed data suggests the adoption of the more parsimonious formulation.
The mathematical formulation adopted for EGP avoids the occurrence of the two problems that emerge with the Steele formulation (Figs 4-7). The typical erroneous negative EGP predictions never happen, due to the intrinsically positive qualitative behaviour of the solutions for EGP. The paradoxical increase in EGP in the first minutes after the start of the experiment also does not occur. Our approach predicts nearly constant production rates after an initial decay immediately after the experimental glucose infusion, with few oscillations and for all subjects, with increasing endogenous glucose productions rates at the end of the experiment, when glucose concentrations tend to their basal values. This behaviour is in fact in excellent accord with the physiological understanding of the processes at play.
Other similar attempts of modelling the endogenous glucose production appeared in the literature [47,48]. In the first work, the Krudys et.al proposed a compartmental model including three differential equations for the EGP production representation (amount of available glucose in the liver, endogenous glucose mass in the accessible compartment and endogenous glucose mass in the slowly equilibrating compartment), with a total of nine parameters plus other two parameters in the equation related to the remote insulin compartment. The model proposed by Krudys et.al therefore appears very expensive from the point of view of parameter estimation, especially in situations where the number of experimental observations is limited. In the work by Visentin et al., EGP is modelled by means of a total of two differential equations and two algebraic equations. The main equation describes the suppression of the EGP as proportional to the glucose rate of change, to glucose concentration and to the delayed insulin action. Two differential equations are used for the description of the dynamics of the delayed insulin action and one algebraic function is used for the contribution to the glucose rate of change, for a total of seven parameters, where one is determined from other model parameters (EGP b ) and two other also appear in other model equations (G b and I b ), for a total of four additional parameters. Besides contemplating an extra parameter to be estimated, it is to be noted that, a-priori, the EGP from Eq (8) can assume negative values due to the fact that the negative terms appearing in Eq (8) are in principle unlimited. This is particularly true for the term proportional to the Glucose derivative which could be very large in correspondence of large glucose administration. A restricted region in the parameter space could be guarantee the positiveness of the solutions and should be investigated.
Our proposed formulation, instead, needs only one additional parameter (λ GI ) to model the glucose contribution by the liver, since the parameter K IG,max can be estimated by steady state conditions. In the present model EGP is represented as an input (an algebraic function) into the plasma glucose compartment, providing a more compact representation of the entire glucose/insulin system.
The estimated values obtained for the insulin sensitivity indices (parameter k xGI,j , j = 1, 2, 3) are of interest. In order to get a better fit of the model to the data, it was necessary to consider three different indices for the three phases of the experimental procedure, each quantifying the effect of insulin on glucose disposal in the corresponding phase. The trend of insulin sensitivity seems to be the same in all subjects: it decreases during the second phase and then increases corresponding to the third OGTT. These results are in agreement with the results of [49], according to which there is intra-day variability of insulin sensitivity. Even if the experiment lasted only 8.30 hours, less than a full day, and even if the variations over time in insulin sensitivity do deserve further ad hoc studies, the results obtained seem to be in line with what was hypothesized in [49]. In [49] the results showed that insulin sensitivity was on average lower at breakfast than at lunch and dinner, despite a great variability between subjects and despite the fact that evidence of a diurnal pattern was demonstrated in a study conducted only on subjects with type I diabetes. It is not possible to say with certainty that the increase observed in this series of experiments is the beginning of a trend leading to greater insulin sensitivity during the last hours of the day. However, the results obtained seem to give credence to this hypothesis.

Conclusion
Although the model appears to fit fairly well with intravenous glucose data, future work will need to address a number of shortcomings. First of all, this approach will need to be adapted to the more complex setting of the estimation of glucose turnover in non-steady state conditions due to the ingestion of glucose. Also, insulin secretion and its effect on glucose disposal should be improved: in fact, while parameters can be reliably estimated under the assumption of autocorrelated error, the reason for such autocorrelation is still obscure. This may imply that the proposed model structure does not represent all necessary physiological mechanisms, or that these are represented in an overly simplistic form. Although improvements can be made, the strength of this model lies in the fact that it respects the physiological assumptions without producing artefacts, as otherwise happens with models which, while producing curves with an apparently better fit to the data, do so at the cost of producing implausible predictions, such as variables that take negative values. In the future it will be therefore necessary to devise a more complex model, which also includes a mathematical representation of the gastrointestinal tract, and because of this it will be necessary to estimate a larger number of parameters, while in any case ensuring that physiology is respected, for instance by exhibiting limited and positive solutions. In conclusion this work shows that appropriate, if simple, modelling of the whole infusion procedure, together with a coherent mathematical representation of the physiology, allows plausible, robust estimation of EGP without recourse to high-frequency sampling or complicated double or triple tracer administration. This approach prevents the artifactual occurrence of negative or spiking EGP predictions, as is instead produced by the Steele procedure.