 Research Article
 Open Access
 Published:
Determination of physiological parameters for endogenous glucose production in individuals using diurnal data
BMC Biomedical Engineering volume 1, Article number: 29 (2019)
Abstract
Background
Triple tracer meal experiments used to investigate organ glucoseinsulin dynamics, such as endogenous glucose production (EGP) of the liver are labor intensive and expensive. A procedure was developed to obtain individual liver related parameters to describe EGP dynamics without the need for tracers.
Results
The development used an existing formula describing the EGP dynamics comprising 4 parameters defined from glucose, insulin and Cpeptide dynamics arising from triple meal studies. The method employs a set of partial differential equations in order to estimate the parameters for EGP dynamics. Tracerderived and simulated data sets were used to develop and test the procedure. The predicted EGP dynamics showed an overall mean R^{2} of 0.91.
Conclusions
In summary, a method was developed for predicting the hepatic EGP dynamics for healthy, prediabetic, and type 2 diabetic individuals without applying tracer experiments.
Background
The plasma glucose concentration in healthy humans is strictly controlled. During fasting conditions the glucose is around 90 mg/dl and after food intake it can reach a concentration of 126144 mg/dl [1]. The glucose comes into the bloodstream by means of absorption from the intestine, breakdown of glycogen (glycogenolysis) and gluconeogenesis [1].
The release of glucose resulting from glycogenolysis and gluconeogenesis is called the endogenous glucose production (EGP). The liver is the major contributor to EGP [2]. Within 12 h after a meal the EGP is suppressed to 30%50% of its fasting value [3–6]. Glucose from the ingested meal will be taken up by the liver, brain, skeletal muscle and adipose tissue in the hours after a meal. The EGP will return to the fasting value after 4 h [7].
Understanding the mechanisms regulating glucose homeostasis is important to better comprehend the effects of new therapies for diseases such as diabetes [8]. Information about EGP dynamics is relevant, because the regulation of this important physiological process is changed in diabetic patients [9–11]. For instance, when fasting hyperglycemia develops in diabetic subjects, elevated rates are observed for both gluconeogenesis and glycogenolysis [12]. Increased gluconeogenesis activity might even be an early feature in the development of glucose dysregulation [13, 14].
Isotope tracer methods are used for the measurement of metabolic (e.g. glucose) flux rates [15], and have been used to investigate amongst others the glucoseinsulin dynamics in organs/tissue (e.g. EGP) [16]. Glucose triple trace isotope labeling experiments can give reliable estimates of EGP, but are expensive and labour intensive [17, 18]. It is therefore not feasible to perform these experiments for large epidemiological and genetic studies. Therefore, we looked for alternative approaches to determine EGP dynamics that avoid the use of tracers.
Several mathematical models describing human glucoseinsulin dynamics have been published in literature [19]. These models are constantly improved or extended to increase their predictive power. Man et al. (20) proposed a glucose dynamics model that describes the glucose and insulin fluxes during a mixed meal. This model is divided into various subsystems, such as gastrointestinal tract, liver, βcell, muscle and adipose tissue. The various parameters of the model in [20] have been estimated from data from glucose isotope labeling experiments in normal humans and humans with type 2 diabetes (T2D). In [20], the EGP in the liver subsystem is described using 2 differential equations and 6 parameters and requires plasma glucose, plasma insulin and portal insulin dynamics data.
Although the model [20] includes parameters settings for normal and type 2 diabetic subject groups, it is of interest to be able to determine these parameters of individual subjects without the need for a new tracer experiment. In this paper, a mathematical procedure is described that is capable of calculating the glucose dynamics model parameters in the model [20] for the liver without glucose isotope labeling experimental data to estimate diurnal EGP dynamics.
Results
The proposed mathematical procedure to estimate the parameters for EGP dynamics was tested on the 5 different data sets. The resulting predicted EGP time profiles were compared to the experimental ones. Residual plots for each data set display the deviation of the predicted EGP.
EGP parameter estimates
Application of the proposed procedure to data from normal, prediabetic and diabetic subjects receiving diurnal mixed meals (data sets 15) resulted in parameter values reported in Table 1.
In particular, the estimated hepatic glucose effectiveness k_{p2} was found to have structurally higher values and the portal insulin signaling parameter k_{p4} structurally lower values compared to the true parameter values for all data sets. The hepatic glucose effectiveness parameter k_{p2} showed similar values when comparing simulated data sets (healthy, both prediabetic conditions (θ_{ds3} and θ_{ds4}) and type 2 diabetic). The value of the estimated parameter k_{p3} in all data sets showed small deviations compared to the optimal parameter values. The parameter A, describing a changed insulin sensitivity during breakfast, showed to be smaller than 1 for all data sets. Parameter differences between the healthy data sets (data sets 1 and 2) and the diabetic data set (data set 5) were observed.
Diurnal EGP dynamics
Figure 1 displays the predictions of diurnal EGP dynamics using the estimated parameter values (θ_{ds15}) (Table 1) vs. the experimental data of data sets 1 to 5. The mean coefficient of determination (R^{2}) for data set 1 was 0.91 and for data set 2 0.91. R^{2} values for prediabetic subjects were 0.97 (data set 3), and 0.87 (data set 4). The R^{2} for a type 2 diabetic subject was 0.92 (data set 5). The less well approximated EGP dynamics were seen in data set 1 and data set 4, consisting of healthy and prediabetic subjects. The EGP dynamics in data set 2 and 3 were most accurately predicted.
Model performance
Figure 2 shows the residual plots for each data set to indicate the model performance. Data set 1 showed the highest error (0.6 mg/kg/min) between the predicted EGP and experimental EGP for breakfast, lunch and dinner. The model showed both under and over estimation for this data set. Lower errors were observed for the data sets 25. The error range of the predicted EGP for data sets (25) was ±0.3 mg/kg/min. Remarkable is that the model systematically under estimated the data from data set 3. On the other hand, the model overestimated the EGP in data set 5.
Reliability of parameter estimates
The parameter values for the data sets 15 were obtained without using the available experimental EGP data, however, this raised the question if these values were estimated reliably compared to the optimal parameter values obtained by fitting the Dalla Man [20] model to the experimental EGP data. To approach this question parameter identifiability analyses employing profile likelihood (PL) estimates were performed. Parameters were considered reliably estimated if they were located within the pointwise confidence interval.
The PL estimates for the parameter values (k_{p2}, k_{p3}, k_{p4} and A) are shown in Fig. 3 for data set 15. The resulting plots reveal practically nonidentifiable and identifiable parameters. However, the observed systematic bias is negligible.
Discussion
This study describes the successful development of a new procedure for the estimation of physiological parameters in a mathematical description of diurnal EGP dynamics in healthy, prediabetic and type 2 diabetic subjects. The procedure exploits the assumption that certain relations between the parameters do not vary between meals during the day. The mathematical description is based on a previously developed computational model describing glucose and insulin dynamics during a mixed meal [20]. A subsystem in this model, duly parameterized, can estimate the EGP using as input the plasma glucose, plasma insulin and portal insulin dynamics data. The various parameters of the subsystem in [20] are estimated from independently determined EGP dynamics, which requires glucose isotope labeling experiments. The proposed procedure obviates the need for such isotope experiments because it can estimate the model parameters directly from data on plasma glucose, insulin and ISR dynamics assembled during 3 consecutive meals.
Several articles reported that the insulin sensitivity is different in the morning (breakfast) compared with the rest of the day or evening [24, 27–29]. To take this diurnal effect into account, we multiplied insulin sensitivity during breakfast vs. the other meals by a factor A. Obtained values of A were between 0.8 and 0.9 for the different data sets. Interestingly, directly fitting the experimental EGP data against the model also yielded an optimal A value smaller than 1 for data set 1.
In the procedure proposed in the present paper, basal EGP value (Eq. 5) was not predicted, but was given to the model as a constant value extracted from the experimental EGP data (i.e., value at t = 0). As a future perspective, a separate prediction procedure could be developed for basal EGP.
The new procedure yielded estimates for the parameters (k_{p2}, k_{p3}, k_{p4} and A) on the basis of which the EGP dynamics were calculated using the model equation [20]. In general the estimated parameter values for the 5 individual data sets showed some deviations compared to the true values. The PL was therefore used to determine the reliability of the obtained parameter values for the different data sets. The analysis demonstrated that a systematic bias is negligible and that most parameters values were reliably estimated because they were located within the defined confidence interval.
One interesting point is that parameter differences between the healthy data sets (data sets 1 and 2) and the diabetic data set (data set 5) were not clearly observed. In fact, while Dalla Man [20] report different values for diabetic vs healthy datasets [20], no standard deviations are reported so one cannot judge the statistical significance of these differences. In the present study, the estimated parameters for data sets 1 and 2 and the diabetic data set (data set 5) may well not differ because they are depending on the particular input data from single individuals. From the PL analysis it can be appreciated that this is indeed the case, as judged from the fairly large confidence intervals around the estimated parameter values. It may still well be possible that when comparing results for groups of individuals, the mean estimated parameter values for the healthy and diabetes groups will show significant differences but this requires a separate, more exhaustive study.
The inaccuracy of the predicted EGP in data set 1, which was extracted from experimental data on healthy subjects, was higher compared to the other data sets (data sets 35), as judged from the residuals and R2 values. The model showed to systematically overestimate the EGP in data set 5. This indicates that the predicted EGP reflects a healthier condition than the experimental EGP. However, the maximal error range is rather small with 0.3 mg/kg/min. Thus, despite the deviations observed for the parameter estimates, the procedure showed a generally accurate prediction of the EGP dynamics for 3 consecutive meals.
In theory, our procedure has broader relevance for a wider spectrum of physiological models. Notably, the procedure might be transferable to models wherein a timedependent variable (here: EGP) cannot be directly measured whereas a deterministic set of equations is available that relates this variable to other timedependent variables (here: G_{p}, I_{d}, and I_{po}) that are measurable.
In this study, five data sets were available for the development of the new and promising methodology to determine EGP dynamics. Only one data set contained experimental diurnal EGP dynamics. However, since the associated experimental errors were not given, it could not be assessed whether the predicted EGP was within the bandwidth set by the experimental inaccuracy. Therefore, it is difficult to precisely judge the quality of the EGP prediction in this procedure. In general, an experimental inaccuracy could be caused by recycling of labeled glucose via glycogenolysis and/or gluconeogenesis, and the equilibration of the tracer/tracee ratio [30]. A next step towards application of the new method in clinical studies will therefore the search for, and possible de novo acquisition, of multiple individual experimental data sets of different subtypes (healthy, prediabetic, T2D subjects), however this is subject of a separate investigation.
Conclusions
In conclusion, a procedure was developed to retrieve physiologically based parameters for the description of diurnal EGP dynamics in subjects without the use of tracer techniques. The procedure makes use of input data consisting of diurnal dynamics of plasma glucose, insulin and ISR. The developed procedure accurately predicted (with an overall mean R^{2} of 0.9) diurnal EGP dynamics data for one experimental and four simulated data sets of healthy, prediabetic and type 2 diabetic subjects. The advantage is that diurnal EGP dynamics can be personally predicted. In the long run, more experimental data from healthy, prediabetic, and type 2 diabetic subjects could contribute to further and extensive validation of the procedure required for clinical applications.
Methods
Computer software
Published data in graphical format was digitized for use in the present paper using the accurate data extraction program GraphClick [21, 22]. All model equations were implemented and analyzed in MATLAB (MATLAB, Version R2012, The MathWorks, Inc., Natick, Massachusetts, United States) (Additional file 2). Ordinary differential equations were simulated using the variable step solver ode15s. The insulin secretion rate (ISR) data were calculated from deconvolution of the Cpeptide time course data using the software program ISEC [23].
Published equations for Endogenous Glucose Production (EGP)
The objective was to develop a method for predicting the parameters of Endogenous Glucose Production (EGP) dynamics in the Dalla Man [20] model from experimental data without needing to apply isotope tracer experiments. The EGP dynamics after a meal in Dalla Man [20] are expressed as a function of glucose and insulin concentrations as schematically shown in Fig. 4.
From this model, the liver subsystem equations were used to formulate a set of equations to calculate EGP dynamics without the need of tracer experiments. From these EGP equations, partial derivatives were formulated and applied to calculate the liver subsystem parameters for estimating diurnal EGP dynamics as explained below.
In the following, the Dalla Man [20] equations (Eqs. 1–5) for the liver subsystem are repeated for convenience. The EGP is described as a function of inputs G_{p} (mg/kg) (glucose mass in plasma), I_{d} (pmol/l) (a delayed insulin signal) and I_{po} (pmol/kg) (insulin in the portal vein) using Eq. 1:
wherein parameter k_{p1} (mg/kg/min) is the extrapolated EGP at zero glucose and insulin, k_{p2} (min^{1}) is liver glucose effectiveness, k_{p3} (mg/kg/min per pmol/l) is a parameter controlling the amplitude of insulin action on the liver, and k_{p4} (mg/kg/min per pmol/kg) is a parameter determining the amplitude of portal insulin action on the liver.
The glucose mass in plasma G_{p} is described using Eq. 2:
wherein V_{g} is the distribution volume of glucose set to 1.88 dl/kg (healthy) or 1.49 dl/kg (type 2 diabetic) and G(t) is the experimental plasma glucose concentration (mg/dl) in time.
The delayed insulin signal, I_{d}, is described using a chain of twocompartments:
wherein I is the plasma insulin concentration (pmol/l) and the rate parameter k_{i} (=0.0079 min^{1} (healthy) or 0.0066 min^{1} (type 2 diabetic) from [20]) determines the delay between the insulin signal and the insulin action.
The portal insulin I_{po}(pmol/kg) dynamics is described by
wherein the pancreatic insulin secretion rate (ISR(t) (pmol/kg/min)) can be calculated from experimental Cpeptide data (details in “Experimental data set: healthy subject” section), and parameter γ (=0.5 min^{1}) is the fixed transfer rate constant between portal vein and liver.
Parameter k_{p1} is calculated as:
where the basal data values for EGP, G_{p}, I_{d}, and I_{po} were extracted from experimental/simulated data sets (details in “Data” section) and k_{p2}, k_{p3}, k_{p4}, and A were estimated as described in the following.
Development of an objective function based on partial derivatives to estimate model parameters
Analyzing the different partial derivatives of Eq. 1 shows interesting relationships. For example, the partial differentiation (\(\frac {\partial {EGP}}{\partial {k_{p2}}}\)) is to first order equal to G_{p}, (\(\frac {\partial {EGP}}{\partial {k_{p3}}}\)) to I_{d}, and (\(\frac {\partial {EGP}}{\partial {k_{p4}}}\)) to I_{po}. This suggests that EGP dynamics might be directly derived from experimental data for G_{p}, I_{d} and I_{po} dynamics. We investigated further based on this interesting observation, by recombining Eq. 1 to find an alternative way to obtain the liver subsystem model parameters i.e. by using experimental data for G_{p}, I_{d}, I_{po} available from diurnal data (3 meals) without the need for isotope tracer experiments.
Taking the partial derivative of Eq. 1, with respect to k_{p2}, k_{p3} and k_{p4}, results in Eqs. 6, 7, and 8, respectively:
where i denotes breakfast (b), lunch (l), or dinner (d).
During the day, biochemical and physiological mechanisms are changing [24]. Especially, a strong variation in glucose tolerance during the day has been observed in different studies [25, 26]. It has also been demonstrated that the insulin responses in healthy subjects are different in the morning (breakfast) compared to the evening (dinner) [27]. Therefore, we decided to introduce a change in the insulin sensitivity during breakfast meal dynamics. To make this assumption more likely to be true, we explicitly took into account the only known intraday variation in EGP regulation, i.e. by introduction of parameter A. Parameter k_{p4} is related to the portal insulin concentration, I_{po}, therefore k_{p4} was multiplied by an additional parameter A only during breakfast dynamics:
where \(k_{p4}^{*} = A \cdot k_{p4}\) only during breakfast meal dynamics (EGP^{b}).
Multiplying Eq. 6 by k_{p2}, Eq. 7 by k_{p3} and Eq. 8 by k_{p4} results in sensitivities, i.e. the change in EGP for normalized changes in the liver parameters, in Eqs. 10, 11, and 12:
Subtracting Eq. 11 from Eqs. 10, 12 from Eqs. 10, and 12 from Eq. 11, results in Eqs. 13, 14 and 15, respectively:
wherein \(C_{23}^{i} = k_{p2} \cdot \frac {\partial {k_{p1}}}{\partial {k_{p2}}}^{i}  k_{p3} \cdot \frac {\partial {k_{p1}}}{\partial {k_{p3}}}^{i}\)
wherein \(C_{24}^{i} = k_{p2} \cdot \frac {\partial {k_{p1}}}{\partial {k_{p2}}}^{i}  k_{p4} \cdot \frac {\partial {k_{p1}}}{\partial {k_{p4}}}^{i}\)
wherein \(C_{34}^{i} = k_{p3} \cdot \frac {\partial {k_{p1}}}{\partial {k_{p3}}}^{i}  k_{p4} \cdot \frac {\partial {k_{p1}}}{\partial {k_{p4}}}^{i}\)
The \(C^{i}_{23}\), \(C^{i}_{24}\), and \(C^{i}_{34}\) values represent mealspecific weighted differences of sensitivities of k_{p1} (extrapolated EGP at zero glucose and insulin) vs. multiple parameter pairs (k_{p2} and k_{p3}, k_{p2} and k_{p4}, and k_{p3} and k_{p4}, respectively).
To estimate single values for k_{p2}, k_{p3}, k_{p4}, and A for a single subject it was assumed that the mathematical constructs called sensitivities (\(C^{i}_{23}\), \(C^{i}_{24}\), and \(C^{i}_{34}\)) do not vary between meals. These mathematical constructs themselves do not translate to simple biological concepts, however they relate parameters that are linked to regulation of EGP. In essence, the assumption is therefore that EGP regulation does not change between 3 successive meals, resulting in:
The model EGP dynamics objective function to obtain the parameter values of k_{p2}, k_{p3}, k_{p4}, and A is therefore formulated by combining Eqs. 16, 17, 18:
wherein i and j are any two combinations of breakfast (b), lunch (l) and dinner (d).
This methodology thus focuses on obtaining estimates only for the parameters which are directly involved in the general EGP dynamics equation (Eq. 1). The parameters k_{i}, V_{g}, and γ relating to glucose kinetics and insulin secretion processes are not directly in the EGP equation and therefore not fitted. However, although these parameters do not appear directly in the equation, they still influence the result of the EGP equation since they multiply/divide the values of measured quantities. An inaccurate assumption of their value could lead to a wrong estimation of the directly involved parameters k_{p2}, k_{p3} and k_{p4}. The model parameter values k_{i}, V_{g}, and γ were extracted from Dalla Man (Table 1) [20]. Here, they obtained these values for normal and type 2 diabetic subjects. So, if the fasting glucose value is below 99 mg/dl (5.5 mmol/l) then the Dalla Man parameter values (k_{i}, V_{g}, and γ) for normal subjects should be used. Parameter values of type 2 diabetic should be used if the fasting plasma glucose value is higher than 126 mg/dl (7.0 mmol/l). For fasting glucose values between 99 mg/dl (5.5 mmol/l) and 126 mg/dl (7 mmol/l) estimates for k_{i}, V_{g}, and γ can be obtained through linear interpolation of Dalla Man [20] normal, and type 2 diabetic subject parameter values.
Briefly, the derivation evolves around a liverrelated parameter k_{p1} which itself explains a rather virtual concept, namely the extrapolated EGP at zero glucose and insulin. This parameter is used to describe the specific behavior of EGP in the virtual condition that no glucose and insulin were present. In the model EGP dynamics objective function, the sensitivity of this parameter to other parameters (i.e. proportional changes of k_{p1} in response to changes in parameters k_{p2}, k_{p3} and k_{p4}) are used. These terms originate in Eqs. 6–8 and appear weighted in Eqs. 10–12. Then, coefficients consisting of differences of pairs of these weighted terms are formed in Eqs. 13–15. The central assumption is that these coefficients do not differ between 3 subsequent meals as expressed in Eqs. 16–18. Hence, the parameter values for a single person are obtained and are equal for breakfast, lunch, and dinner. As stated earlier, the essence of the underlying biological hypothesis is that the parameters for regulation of hepatic EGP production do not vary during 3 consecutive meals within a single day, but the derivation follows a mathematical rather than a biological reasoning.
Data
The procedure was applied to an experimental data set and simulated data sets (5 data sets in total) representing different health conditions of subjects as follows:
Experimental data set: healthy subject
Data set 1: Healthy subject Experimental mean 3meal data curves of 20 healthy subjects (i.e. mean plasma glucose, plasma insulin, Cpeptide and EGP) were extracted from Fig 1A, 1B, and 2A in [28] using GraphClick [22]. The data were converted to the desired units (Glucose: mmol/l to mg/dl (×18.02); EGP: μmol/kg/min to mg/kg/min (×0.18)). Each segment (breakfast, lunch and dinner), starts when the meal is given at time zero, and ends 240 minutes later. The ISR was calculated from Cpeptide data (Fig 1B in [28]) using the software program ISEC [23].
Preprocessing of the data set was necessary in order to obtain a single set of plasma glucose, plasma insulin and ISR data curves for a hypothetical representative individual subject in the study of [28] that were consistent with the mean EGP curves of the 20 subjects. We adopted this procedure because the individual data sets remained unavailable. The procedure is explained and the raw data can be seen in (Additional file 1). The result, i.e. the final input data for the model, is depicted in Fig. 5. Data was collected at time points 0, 15, 30, 45, 60, 90, 120, 150, 180, 210, 240 minutes. This data set is referred to as data set 1 throughout the paper.
Simulated data sets: healthy, prediabetic, and diabetic subjects
Data set 2: Healthy subject Simulated 3meal data curves (i.e. plasma glucose, plasma insulin, ISR and EGP) for a healthy subject were generated using the meal simulation model in [20]. To this end, the in silico model in [20] was implemented in MATLAB and the parametric portrait of a normal subject (Table 1 in [20]) was used for simulation.The predicted plasma glucose, plasma insulin, ISR and EGP were based on 3 consecutive meals i.e. a breakfast (8 a.m.) meal of 45 gram, a lunch (12 p.m.) meal of 70 gram, and a dinner (8 p.m.) meal of 70 gram. Each meal segment starts when the meal is given (t=0 for that segment) and ends 240 min later. Data was collected at time points 0, 15, 30, 45, 60, 90, 120, 150, 180, 210, 240 min. The data is referred to as data set 2 throughout the paper, and can be seen in Fig. 6.
Data sets 3 and 4: Prediabetic subjects The in silico model in [20] implemented in MATLAB was used to obtain simulation data for 2 prediabetic subjects having different subtypes i.e.(1) glucose intolerance, and (2) decreased insulin action. For the glucose intolerant subject, the parameter values for normal subjects ([20], Table 1) were used, except that V_{max} and k_{p3}, were halved. For the subject with decreased insulin action, parameter values for normal subjects were used except that K and β were halved. The predicted plasma glucose, plasma insulin, ISR and EGP were based on 3 consecutive meals i.e. a breakfast (8 a.m.) meal of 45 gram, a lunch (12 p.m.) meal of 70 gram, and a dinner (8 p.m.) meal of 70 gram. Each meal segment starts when the meal is given (t=0 for that segment) and ends 240 min later. Data was collected at time points 0, 15, 30, 45, 60, 90, 120, 150, 180, 210, 240 min. The data are referred to as data set 3 (glucose intolerant subject) and data set 4 (subject with a decreased insulin action) throughout the paper, and can be seen in Figs. 7 and 8.
Data set 5: type 2 diabetic subject The in silico model in [20] implemented in MATLAB was used to obtain simulation data for a type 2 diabetic subject. The parameter values for type 2 diabetic subjects ([20], Table 1) were used. The predicted plasma glucose, plasma insulin, ISR and EGP were based on 3 consecutive meals i.e. a breakfast (8 a.m.) meal of 45 gram, a lunch (12 p.m.) meal of 70 gram, and a dinner (8 p.m.) meal of 70 gram. Each meal segment starts when the meal is given (t=0 for that segment) and ends 240 min later. Data was collected at time points 0, 15, 30, 45, 60, 90, 120, 150, 180, 210, 240 min. The data are referred to as data set 5 throughout the paper, and can be seen in Fig. 9.
Parameter optimization strategy
The EGP data fitting procedure was implemented in Matlab with the function lsqnonlin, minimizing the model EGP dynamics objective function (Eq. 19). Lower boundaries were set to 0.0009, 0.001, 0.01, 0.8 and upper boundaries to 0.02, 0.05, 0.09, 1.4 for k_{p2}, k_{p3}, k_{p4} and A, respectively. The function multistart was used to start the lsqnonlin function from multiple start points (150) in order to increase the probability of finding a global optimum. The reported parameters for each data set were those for which the model EGP dynamics objective function achieved its minimal value.
The obtained parameter values for each data set were also checked against the optimal values, which are determined from directly fitting the experimental EGP against the model. These optimal parameter values were calculated by minimizing the following cost function (V) with the above explained settings:
wherein the cost function V is defined as the squared difference between the model output (\(EGP^{calc}_{i,j}\)) and experimental EGP (\(EGP^{exp}_{i,j}\)) for every time step, i, the model parameters p, and for breakfast, lunch and dinner, j. Here 5 data sets that originated from the experimental and simulated data (Figs. 5, 6, 7, 8 and 9) were treated as experimental data (\(EGP^{exp}_{i,j}\)).
Parameter identifiability analysis
To investigate the reliability of the obtained parameter estimates (k_{p2}, k_{p3}, k_{p4} and A), parameter identifiability analysis using PL [31] was performed. PL allows to derive the identifiability of parameters in nonlinear models, to design optimal experiments that improve parameter identification and to calculate likelihoodbased confidence intervals.
First, the weighted sum of squared residuals was set as the PL objective function measuring the agreement of experimental data with modelpredicted data:
where \(y^{exp}_{kl}\) denotes the experimental EGP data for each observable time point k, \(\sigma ^{exp}_{kl}\) is the measurement error of the experimental data and was assumed to be normally distributed. Here 5 data sets that originated from the experimental and simulated data (Figs. 5, 6, 7, 8 and 9) were treated as experimental data \(\left (y^{exp}_{kl}\right)\).
The \(y^{cal}_{k}\) are calculated EGP values at time points when experimental data were measured, m is the number of time points and d is the number of meal segments.
In the PL procedure, the parameters were calculated numerically:
where the value \(\chi ^{2}(\hat {\theta })\) corresponds to the maximum likelihood estimation for Gaussian measurement noise. The PL is determined by an iterative procedure where the first of the parameters is slightly shifted from its optimal value [31, 32] whereafter the other parameters are fitted again. Then the parameter is shifted again and the other parameters are fitted again. The shifting is done both towards smaller values, and towards higher values. This procedure for a given parameter will end in either direction when a threshold in the likelihood is met. After that, the procedure is repeated for the other parameters. The likelihoodbased confidence interval for a parameter was calculated to define a confidence region [31, 33]:
with df the number of degrees of freedom and α the confidence level. The df gives for df=1 the pointwise and for df=4 (number of fitted parameters) the simultaneous confidence interval, both to the confidence level α of 0.68, which corresponds to one standard deviation.
A parameter is identifiable if the confidence interval [ σ^{−}, σ^{+}] of its estimate \(\hat {\theta }\) is finite. This indicates that this parameter can be calculated from experimental data sets. Practical nonidentifiability occurs when a parameter shows a likelihoodbased confidence region which is infinite extended in the direction of θ, even if the likelihood has a unique minimum for this parameter [31]. The confidence interval of a practically nonidentifiable parameter is not automatically infinitely extended to both sides; a finite lower or upper bound of the confidence interval [ σ^{−}, σ^{+}] can still exist. Resolving a practically nonidentifiable parameter is possible by increasing the amount and quality of experimental data and/or the choice of the time points [31]. Structural nonidentifiability occurs if the confidence interval is infinite to both sides. Structural nonidentifiability cannot be solved by a more accurate experimental data set, and indicates that there is missing information [31].
With this method, obtained parameter values were inspected on the deviations from \(\chi ^{2}(\hat {\theta })\) and if the corresponding predicted model observable (\(y^{cal}_{kl}\)) was in agreement with experimental noise, assumed to be \(y^{exp}_{kl}*0.1 + \text {max}(y^{exp}_{kl})*0.05\), so as to gain information on the quality of the model predictions.
Availability of data and materials
All data and software generated or analysed during this study are included in this published article and its supplementary information files.
Abbreviations
 df:

Degrees of freedom
 EGP:

Endogenous glucose production
 ISR:

Insulin secretion rate
 PL:

Profile likelihood
 T2D:

Type 2 diabetes
References
 1
Frayn KN. Insulin resistance and lipid metabolism. Curr Opin Lipidol. 1993; 4(3):197–204.
 2
Gerich J. Role of the kidney in normal glucose homeostasis and in the hyperglycaemia of diabetes mellitus: therapeutic implications. Diabet Med. 2010; 27(2):136–42.
 3
Frayn KN. Metabolic Regulation: A Human Perspective. Chichester: Wiley; 2009.
 4
Capaldo B, Gastaldelli A, Antoniello S, Auletta M, Pardo F, Ciociaro D, et al.Splanchnic and leg substrate exchange after ingestion of a natural mixed meal in humans. Diabetes. 1999; 48(5):958–66.
 5
Singhal P, Caumo A, Carey PE, Cobelli C, Taylor R. Regulation of endogenous glucose production after a mixed meal in type 2 diabetes. Am J PhysiolEndocrinol Metab. 2002; 283(2):E275–83.
 6
Woerle HJ, Meyer C, Dostou JM, Gosmanov NR, Islam N, Popa E, et al.Pathways for glucose disposal after meal ingestion in humans. Am J PhysiolEndocrinol Metab. 2003; 284(4):E716–25.
 7
Fery F, d’Attellis N, Balasse E. Mechanisms of starvation diabetes: a study with double tracer and indirect calorimetry. Am J PhysiolEndocrinol Metab. 1990; 259(6):E770–7.
 8
Beysen C, Hellerstein MK, Turner SM. Isotopic Tracers for the Measurement of Metabolic Flux Rates. In: Translational Research Methods for Diabetes, Obesity and Cardiometabolic Drug Development. Springer: 2015. p. 71–97. https://doi.org/10.1007/9781447149200_3.
 9
Bogardus C, Lillioja S, Howard B, Reaven G, Mott D. Relationships between insulin secretion, insulin action, and fasting plasma glucose concentration in nondiabetic and noninsulindependent diabetic subjects. J Clin Investig. 1984; 74(4):1238.
 10
Rizza RA. Pathogenesis of fasting and postprandial hyperglycemia in type 2 diabetes: implications for therapy. Diabetes. 2010; 59(11):2697–707.
 11
Perriello G, Pampanelli S, Del Sindaco P, Lalli C, Ciofetta M, Volpi E, et al.Evidence of increased systemic glucose production and gluconeogenesis in an early stage of NIDDM. Diabetes. 1997; 46(6):1010–6.
 12
Gastaldelli A, Baldi S, Pettiti M, Toschi E, Camastra S, Natali A, et al.Influence of obesity and type 2 diabetes on gluconeogenesis and glucose output in humans: a quantitative study. Diabetes. 2000; 49(8):1367–73.
 13
Bock G, Chittilapilly E, Basu R, Toffolo G, Cobelli C, Chandramouli V, et al.Contribution of hepatic and extrahepatic insulin resistance to the pathogenesis of impaired fasting glucose role of increased rates of gluconeogenesis. Diabetes. 2007; 56(6):1703–11.
 14
Basu R, Barosa C, Jones J, Dube S, Carter R, Basu A, et al.Pathogenesis of prediabetes: role of the liver in isolated fasting hyperglycemia and combined fasting and postprandial hyperglycemia. J Clin Endocrinol Metab. 2013; 98(3):E409–17.
 15
Wolfe RR, Chinkes DL. Isotope tracers in metabolic research: principles and practice of kinetic analysis. Chichester: Wiley; 2005.
 16
Vella A, Rizza RA. Application of isotopic techniques using constant specific activity or enrichment to the study of carbohydrate metabolism. Diabetes. 2009; 58(10):2168–74.
 17
Dalla Man C, Caumo A, Basu R, Rizza R, Toffolo G, Cobelli C. Minimal model estimation of glucose absorption and insulin sensitivity from oral test: validation with a tracer method. Am J Physiol Endocrinol Metab. 2004; 287(4):E637–43.
 18
Rizza RA, Toffolo G, Cobelli C. Accurate Measurement of Postprandial Glucose Turnover: Why Is It Difficult and How Can It Be Done (Relatively) Simply?Diabetes. 2016; 65(5):1133–45.
 19
Ajmera I, Swat M, Laibe C, Le Novere N, Chelliah V. The impact of mathematical modeling on the understanding of diabetes and related complications. CPT: Pharmacomet Syst Pharmacol. 2013; 2(e54):1–14.
 20
Man CD, Rizza RA, Cobelli C. Meal Simulation Model of the GlucoseInsulin System. Biomed Eng IEEE Trans. 2007; 54(10):1740–9. https://doi.org/doi:10.1109/TBME.2007.893506.
 21
Boyle MA, Samaha AL, Rodewald AM, Hoffmann AN. Evaluation of the reliability and validity of GraphClick as a data extraction program. Comput Hum Behav. 2013; 29(3):1023–7. http://dx.doi.org/10.1016/j.chb.2012.07.031.
 22
GraphClick (Version 3.0.2). 2010. http://www.arizonasoftware.ch. Accessed 08 Jan 2018.
 23
Hovorka R, Soons PA, Young MA. ISEC: a program to calculate insulin secretion. Comput Methods Prog Biomed. 1996; 50(3):253–64.
 24
Maury E, Ramsey KM, Bass J. Circadian rhythms and metabolic syndrome from experimental genetics to human disease. Circ Res. 2010; 106(3):447–62.
 25
Carroll KF, Nestel PJ. Diurnal variation in glucose tolerance and in insulin secretion in man. Diabetes. 1973; 22(5):333–48.
 26
Van Cauter E, Blackman JD, Roland D, Spire JP, Refetoff S, Polonsky KS. Modulation of glucose regulation and insulin secretion by circadian rhythmicity and sleep. J Clin Investig. 1991; 88(3):934.
 27
Lee A, Ader M, Bray GA, Bergman RN. Diurnal variation in glucose tolerance: cyclic suppression of insulin action and insulin secretion in normalweight, but not obese, subjects. Diabetes. 1992; 41(6):750–9.
 28
Saad A, Dalla Man C, Nandy DK, Levine JA, Bharucha AE, Rizza RA, et al.Diurnal pattern to insulin secretion and insulin action in healthy individuals. Diabetes. 2012; 61(11):2691–700.
 29
Boden G, Ruiz J, Urbain JL, Chen X. Evidence for a circadian rhythm of insulin secretion. Am J Physiol Endocrinol Metab. 1996; 271(2):E246–52.
 30
Tigas SK, Sunehag AL, Haymond MW. Impact of duration of infusion and choice of isotope label on isotope recycling in glucose homeostasis. Diabetes. 2002; 51(11):3170–5.
 31
Raue A, Kreutz C, Maiwald T, Bachmann J, Schilling M, Klingmüller U, et al.Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Bioinformatics. 2009; 25(15):1923–9.
 32
Venzon D, Moolgavkar S. A method for computing profilelikelihoodbased confidence intervals. Appl Stat. 1988:87–94. https://doi.org/10.2307/2347496.
 33
Meeker WQ, Escobar LA. Teaching about approximate confidence regions based on maximum likelihood estimation. Am Stat. 1995; 49(1):48–53.
Acknowledgments
Not applicable.
Funding
This study was partially funded by the European Commission under the 7th Framework Programme (MISSIONT2D project, contract No.600803). The funding had no role in the design of the study and collection, analysis, and interpretation of data and writing the manuscript.
Author information
Affiliations
Contributions
MFvS and SK conceived the study and performed the analysis. MFvS, AAdG, AKG, and SK analyzed the results. AAdG, AKG and MFvS prepared the manuscript. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Preprocessing of data set 1 (experimental data set).
Ref. [28] gives mean data input curves of 20 subjects. To test the new procedure developed in this paper, a set of data curves for plasma glucose ( G_{p}), plasma insulin ( I_{p}) and insulin secretion rate (ISR) for a single subject needed to be generated that would produce the mean EGP data curves of the 20 subjects. To this end, the parameters (p), k_{p2}, k_{p3} and k_{p4} plus additional parameters representing the mean input data (i.e. G_{p}, I_{p} and ISR) at time points 20, 30, 60, and 90 minutes for 3 meals (see below), were optimized using the gradientbased least squares solver lsqnonlin in MATLAB (see Additional file 1: Figure S1). This was performed by minimizing the preprocessing experimental data objective function (Eq. 24):
\( Obj (p) = \sqrt {(EGP_{i,j}^{calc}(p,t)  EGP_{i,j}^{exp}(t))^{2}}, \qquad \text {(24)} \)
wherein the function Obj is defined as the squared difference between the Dalla Man [20] model output (\(EGP^{calc}_{i,j}\)) and experimental EGP (\(EGP^{exp}_{i,j}\)) for every time step, i, and for breakfast, lunch and dinner, j. For experimental data (\(EGP^{exp}_{i,j}\)) in Eq. 24 is taken from Saad et al. [28].
The additional parameters for singlesubject G_{p}, I_{p}, and ISR at t=20, 30, 60, and 90 minutes were defined as the shifts from the mean data values and were constrained within ±10 mg/dl, ±25 pmol/l, and ±1 pmol/kg/min, respectively. Additional file 1: Table S1 displays the values for these shifts that resulted from the optimization procedure. The resulting singlesubject data were used as data set 1 (Fig. 2).
Additional file 2
mATLAB code. MATLAB code of the described methodology.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
van Stee, M.F., Krishnan, S., Groen, A.K. et al. Determination of physiological parameters for endogenous glucose production in individuals using diurnal data. BMC biomed eng 1, 29 (2019). https://doi.org/10.1186/s424900190030z
Received:
Accepted:
Published:
Keywords
 Endogenous glucose production (EGP)
 Nonlinear system
 Parameters
 Diurnal information
 Partial differentiation