- Email: [email protected]
ANALYSIS OF DOUBLE-TRACER GLUCOSE KINETICS IN HUMANS DURING ORAL GLUCOSE TOLERANCE TEST
ANALYSIS OF DOUBLE-TRACER GLUCOSE KINETICS IN HUMANS DURING ORAL GLUCOSE TOLERANCE TEST Karl Thomaseth ∗ Amalia Gastaldelli ∗∗ Alessandra Pavan ∗ Rachele Berria ∗∗∗ Leonard Glass ∗∗∗ Ralph DeFronzo ∗∗∗ ∗
CNR Institute of Biomedical Engineering, Padova, Italy ∗∗ CNR Institute of Clinical Physiology, Pisa, Italy ∗∗∗ Univ of Texas Health Science Center at San Antonio, USA
Abstract: A mathematical model based on minimal physiological assumptions for describing simultaneously multiple glucose measurements during modified oral glucose tolerance test is presented. Parameter identification was carried out using a population approach, which allowed precise characterisation of average glucose kinetic parameters in the studied cohort, as well as between–subject variability associated with glucose tolerance state and other covariates. The statistical and computational complexity added with the formulation of non-linear hierarchical population kinetic models is compensated by improved robustness, over the traditional two–stage approach, in selecting the most adequate model structure/order at the individual level and the most significant determinants of c 2006 IFAC glucose kinetics at the population level. Copyright Keywords: Mathematical models, Medical applications, Parameter identification, Physiological models, Reduced-order models, Statistical analysis
1. INTRODUCTION The oral glucose tolerance test (OGTT) is increasingly used as clinical test for studying glucose tolerance and insulin secretion in cohorts at risk of type 2 diabetes and for assessing the efficacy of new antidiabetic drugs. In its simplest form (oral load of 75 g glucose dissolved in water and a single measure of plasma glucose concentration after 2 hours) the OGTT is recommended by the World Health Organization for the diagnosis of diabetes in presence of elevated fasting plasma glucose. Given its moderate invasiveness for patients and usefulness in providing dynamic information on the kinetics of glucose and other endogenous substances, more complex OGTT-based protocols have been employed, e.g. with frequent blood sampling and administration of exogenous sub-
stances, combined with mathematical modelling of dynamic concentration profiles. In this study we consider the double-tracer labeled OGTT, which consists of a primedcontinuous infusion of tritiated glucose ([3-3 H]glucose) (to quantify endogenous glucose turnover by achieving steady state concentrations of [33 H]-glucose) and carbon labeled glucose ([1-14 C]glucose) added to the oral glucose load (to determine exogenous glucose absorption rates). The applicability of mathematical modelling to these kinds of experiments have been already demonstrated by others (Dalla Man et al., 2002; Dalla Man et al., 2004; Dalla Man et al., 2005). In particular, the so-called Bergman’s minimal model of glucose disappearance (MINMOD) (Bergman et al., 1979), which had been developed for modelling
frequently sampled IVGTT (FSIGT) data, has proven to adequately describe, after straightforward modifications, also glucose kinetics during OGTT. More specifically, by explicitly modelling the exogenous glucose appearance rate (Ra ) of orally administered unlabeled and [1-14 C]-glucose, it has been shown that MINMOD is able to accurately match the Ra profile determined with another, independent, modelling approach and to provide a figure of insulin sensitivity of glucose disappearance equivalent to that measured with the more laborious glucose clamp technique. Disadvantages of the model proposed in (Dalla Man et al., 2005) include the relatively high number of model parameters to be estimated from experimental data collected in a single individual and the use of two different models (similar structure but with different parameter values) for unlabeled and tracer glucose kinetics. Aims of this study were to re-assess the minimal modelling approach applied to double-tracer OGTT experiments within a unifying framework using: (i) a single model for describing simultaneously [3-3 H]-glucose, [1-14 C]-glucose and unlabeled glucose kinetics; (ii) a simplified model with reduced number of parameters, each having a unique value and physiological interpretation for different glucose measurements; (iii) a more robust population parameter estimation approach to better separate statistical variability of parameter estimates due to measurement and modelling errors from between-subject variability of kinetic parameters dependent from different metabolic states and other covariates, the characterisation of which is of major interest in population studies. The next section briefly summarises the studied cohort and the experimental protocol, and provides a detailed description of the mathematical model used to fit experimental data and the statistical population approach used to estimate model parameters. Section 3 presents results of the population modelling analysis, which are commented in section 4 together with a critical discussion about the proposed modelling approach.
2. MATERIALS AND METHODS 2.1 Subjects and protocol Twelve type 2 diabetic patients (7 males and 5 females; age = 53.6±2.5 years; body weight (BW) = 81.9±3.6 kg; body mass index (BMI) = 30.5±1.1 kg m−2 ) and ten normal glucose-tolerant subjects (5 males and 5 females; age = 39.6±3.7 years; BW = 90.1±3.9 kg; BMI = 31.1±0.9 kg m−2 ), matched for BMI and used as control subjects, underwent after overnight fast a doubletracer OGTT consisting of a primed-constant in-
fusion of [3-3 H]-glucose (0.25 µCi min−1 ) starting at time −120 min and an oral glucose load at time 0 of 75 g of glucose diluted in water containing 75 µCi [1-14 C]-glucose. Blood samples were collected every 15 minutes starting at time −30 min until 240 min for determination of plasma concentrations of unlabeled and labeled glucose and insulin. Concentration measurements are expressed as follows: glucose (mg dl−1 ), [3-3 H]glucose (dpm ml−1 ), [1-14 C]-glucose (dpm ml−1 ) and plasma insulin (µU ml−1 ).
2.2 Mathematical Modelling The model of glucose kinetics was formulated by assuming, for any of the three different glucose measurements, a one-compartment structure with fixed unknown glucose distribution volume (VG ). The generic mass balance equation is described by ˙ X(t) = −Rd (t) + Rp (t) + Ra (t)
(1)
where X(t) represents compartmental mass, expressed in (mg) and (dpm) for unlabeled and ˙ tracer glucose, and with time derivative X(t) (mg min−1 ) and (dpm min−1 ), respectively; glucose disappearance rate (Rd ) is defined as the product of instantaneous glucose clearance (CL) (dl min−1 ) and the corresponding plasma glucose concentration; the endogenous glucose production rate (Rp ) is nonzero only for unlabeled glucose; the exogenous glucose appearance rate (Ra ) is constant and known for [3-3 H]-glucose and timevarying and unknown for unlabeled and [1-14 C]glucose, yet with known total area under the curve (AUC) over an infinite horizon equal to the administered doses. Modelling glucose disappearance rate (Rd ). Given the above definition of Rd it is sufficient to define the mathematical structure for glucose clearance (CL). According to the minimal model of glucose disappearance (MINMOD) CL is described as the sum of a constant insulin-independent term, the glucose effectiveness at zero insulin (GEZI) (Abbate et al., 1993), and a time-varying term given as the product of insulin sensitivity index (SI ) (dl min−1 µU−1 ml) and insulin concentration in a remote compartment (IR (t)) CL(t) = GEZI + SI IR (t)
(2)
It is worth noting that the adopted parameterisation is different, but mathematically equivalent, to that conventionally used in the MINMOD literature, i.e. parameters expressing clearances (dl min−1 ) are obtained by multiplying the conventional counterpart of fractional clearances (min−1 ) by the glucose distribution volume (VG );
moreover the more commonly used glucose effectiveness at basal insulin (SG ) is related to the other model parameters through the equation SG = GEZI + SI Ib
(3)
The dynamics of insulin concentration in a remote compartment is described by I˙ R (t) = −P2 (IR (t) − I(t))
(4)
tion rate (Rp ) being dependent on deviations of remote insulin from basal through an empirical Hill-type equation. Given the glucose effectiveness at basal insulin (SG ), Eq. (3), and basal glucose concentration (Gb ), endogenous glucose production rate was modeled as Rp (t) = SG Gb
I50 h I50 h + (IR (t) − Ib )h
(8)
with initial condition equal to basal plasma insulin concentration (Ib ); plasma insulin concentration profile (I(t)) reconstructed by linear interpolation of measurements; and remote insulin turnover rate (P2 ). Thus IR (t) is simply a lagged profile of plasma insulin, and is related to the traditional MINMOD insulin action by: X(t) = SI (IR (t) − Ib ).
where I50 characterises the remote insulin concentration for 50% inhibition of basal glucose production, and the exponent h is an empirical parameter.
Modelling exogenous glucose appearance rate (Ra ). Exogenous unlabeled and [1-14 C]-glucose, administered through the oral load, share the same absorption time profile. A mathematical function, Ra (t, p), parametrised by p, is a suitable candidate for representing gastrointestinal glucose absorption if it is nonnegative and its total AUC is normalised to one, such that multiplication by administered dose yields effective absorption rate. For this purpose we adopted the following definition for Ra (t, p):
G ˙ G (t) = −CL(t) X (t) + RG X p (t) + Ra (t) (9) VG
T
T
β
Ra (t, p) dt = 1 − e−( α )
(5)
0
which ensures, if the two parameters α and β are positive, that the incremental AUC of glucose absorption strictly increases with time and tends exponentially to 1 for T → ∞. Given the oral unlabeled glucose dose (DG ) and the specific activity of administered [1-14 C]14 glucose (Sa C ), the exogenous glucose appearance rate for unlabeled and [1-14 C]-glucose are: G RG a (t) = D Ra (t, p) 14
C
14
Ra (t) = Sa
C
G
D Ra (t, p)
(6) (7)
Modelling endogenous glucose production rate (Rp ). In MINMOD endogenous glucose production is not modeled explicitly, but it is assumed that the control of net hepatic glucose balance (NHGB) parallels that of peripheral glucose utilisation. The MINMOD-derived insulin sensitivity index (SI ) for unlabeled glucose accounts thus also for the effects of insulin on NHGB. This leads to different values if SI is estimated from tracer glucose kinetics, which depends on glucose disposal alone. In view of our aim of defining a unique sensitivity index SI we assumed endogenous glucose produc-
Final model equations. The mathematical model equation used to fit double-tracer OGTT experiments can be summarised as follows G
14
˙ 14 C (t) = −CL(t) X X
14 (t) + Ra C (t) VG
3
˙ 3 H (t) = −CL(t) X X
C
H
(t)
VG
(10)
3
H + uinf
(11) 14
with initial conditions XG (0) = Gb VG , X C (0) = 3 3 0, X H (0) = GbH VG where CL(t) is derived from Eq. (2) with remote insulin determined through G Eq. (4); RG p (t) is expressed by Eq. (8); Ra (t) 14 3 C H and Ra (t) by Eq. (6) and (7), respectively; uinf represents the steady state constant infusion rate. The final set of estimated model parameters consist of: {P2 , GEZI, SI , VG , α, β, I50 , h}, with the first four parameters pertaining to the original MINMOD and two new parameters each for describing glucose absorption and production. Additional assigned model parameters are: test dose 14 3 H DG (75 g); individually measured Sa C and uinf ; 3 H basal concentration values Gb , Ib and Gb obtained for each individual as the average of measurements taken before the beginning of the test. Moreover, individual plasma insulin concentration profile (I(t)) is determined from insulin concentration measurements. Model equations were implemented using the modelling software Pansym (Thomaseth, 2003).
2.3 Parameter estimation Having defined the kinetic model for describing the concentration profile of unlabeled glucose, [33 H]-glucose and [1-14 C]-glucose, the i-th data set of the i-th subject can be formally represented through a non-linear regression model such as yi = f (ti , θi ) + ei
(12)
where yi is the vector of all measurements collected in the i-th subject, f (ti , θi ) represents the model predictions at the sampling times ti for given vector of parameters θi , and ei represents measurement (and modelling) noise, which was chosen among various alternatives as homoscedastic with a different variance estimated for each type of glucose measure. To guarantee the fulfillment of positivity constraints on parameters and to reduce the risk of identification problems in presence of large interindividual variability of kinetic parameters, logarithmic transformations were applied to model parameters before their estimation. Thus θ represent in the following the vector of log-transformed parameters, i.e. {log P2 , log GEZI, log SI , log VG , log α, log β, log I50 , log h}, which were the actually estimated parameters. Before their use in numerical simulation of model equations they were backtransformed through exponentiation. Results and statistical figures of precision (confidence intervals) will also be presented for the original, backtransformed, parameters. For obtaining information on the average kinetic response and on between-subject variability, indiˆ i , could have been vidual parameter estimates, θ obtained for each subject by non-linear weighted least squares fitting with subsequent statistical analysis to determine those covariates that affect kinetic parameters (this is the so-called twostep procedure). A more robust and dependable approach is the population modelling technique based on simultaneous analysis of all data using a statistical model that explicitly accounts for intra- and inter-individual variations of parameters (Davidian and Giltinan, 1995). In particular the non-linear mixed effects (NLME) approach was employed, which is based on maximum likelihood estimation of model parameters and covariance matrices of random variables under the assumption of normal distributions (Pinheiro and Bates, 2000). In this context, the k-th component of the individual parameter vector θi is in general described by a linear model as follows θik = θ0k + Xik ξ k + δik
(13)
where θ0k represents the average population value, Xik is a row vector of explicative variables (continuous and factorial covariates) to predict interindividual variability. Parameters θ0k and ξ k are called fixed effects because, by hypothesis, they don’t change between subjects of the studied population; δik are the so-called random effects that represent unpredictable between-subject variations. Yet, not all kinetic parameters were described using the complete structure of Eq. (13). To select the minimum number of estimated pa-
rameters able to quantify the statistically significant determinants of intra- and inter-individual variability, the inclusion of a covariate within Xik or the actual need for a random effect δik was assessed on the basis of physiological plausibility of results and statistical criteria for precision of parameter estimates and minimum model complexity (Davidian and Giltinan, 1995). For instance, the decision rule for including a covariate in the model was that the corresponding multiplicative regression parameter be estimated with a significance level of at least P=0.05 of being different from zero.
3. RESULTS The ability of the model to fit experimental data was excellent on pooled data (not shown) and very good on average data for normal control and diabetic groups, separately. Figure 1 shows the average concentrations (mean±SE) and model predictions for cold glucose, [3-3 H]-glucose, and [114 C]-glucose. For direct comparison each matching glucose measure is plotted with same vertical range but possibly with different offset. With regard to model performance at an individual level, the percentage error in model prediction (residual/fitted) had a low bias, being on average -0.02% for unlabeled glucose, 2.1% for [3-3 H]glucose, and 0.47% for [1-14 C]-glucose, and had an acceptable dispersion with standard deviations (a measure of mean coefficient of variation) of 6.2%, 16.4% and 7.8%, respectively. Fixed effects parameter estimates, characterizing both subgroups at the population level as one ensemble, are reported in Table 1 together with figures of precision, i.e. coefficient of variation and lower (L) and upper (U) limits of 95% statistical confidence intervals. It can be observed that parameter estimates have physiologically plausible values, such as GEZI = 0.72 dl min−1 that corresponds to a fractional glucose clearance of 0.0067 min−1 ; SI = 0.098 dl min−1 µU−1 ml in normal subjects, corresponds to a value of 9.17 10−4 min−1 µU−1 ml. The insulin sensitivity index in diabetic patients resulted in low values (only 4.9% of SI of normal subjects). In addition to the differences in SI between normal and diabetic subjects the only covariate found to significantly affect glucose kinetics was BW influencing GEZI. The value of %GEZI/kg reported in Table 1 quantifies the % variation in glucose effectiveness at zero insulin for increases in 1 kg BW. The reported reference value of GEZI corresponds to an arbitrarily chosen nominal BW of 85 kg. An interesting result is the quite low value estimated of 16 µU ml−1 for insulin concentration for 50% inhibition of basal glucose production (I50 ).
100
150
200
350 250
Glucose (mg/dl) 50
150
150
Glucose (mg/dl)
0 50 0
250
0
50
50
100
150
200
250
150
200
250
time(min)
200
250
200
250
6500
50
100
150
5000 2000 0
[1-14C]-Glucose (dpm/ml)
5000 2000
100
250
time(min)
0
[1-14C]-Glucose (dpm/ml)
50
200
5500 0
time(min)
0
150
4500
[3-3H]-Glucose (dpm/ml)
3500 0
100
time(min)
2500 1500
[3-3H]-Glucose (dpm/ml)
time(min)
0
50
100
150
time(min)
Fig. 1. Experimental data (mean±SE) and best average model predictions in normal controls (left panels) and diabetic subjects (right panels). Equivalent data are drawn with same scales but, whenever convenient, with different offsets. Table 1. Fixed effects population parameter estimates and precision. Parameter P2 GEZI %GEZI/kg SInormal SIdiabetic VG α β I50 h
Value 0.0163 0.72 2.09 0.0981 0.0048 107 123.1 1.451 16.11 1.425
CV% 28.9 9.74 0.70 14.8 18.2 18.1 9.3 3.1 24.7 22.7
95%(L) 0.0093 0.59 0.71 0.073 0.0033 75.1 102.5 1.36 9.92 0.91
95%(U) 0.0287 0.87 3.5 0.13 0.0068 152.5 147.8 1.55 26.2 2.22
4. DISCUSSION This is the first report, to our knowledge, on simultaneous modelling of multiple sampled glucose data from double-tracer oral glucose tolerance test in healthy and diabetic subjects using a population modelling approach. A similar approach has been adopted for the analysis of nonsteady state glucose kinetics during euglycemic hyperinsulinemic glucose clamp performed in cats (Hoenig et al., in press). The proposed model is closely related to previous work by others, e.g. (Dalla Man et al., 2004), the main difference being that glucose production is modeled here explicitly with a heuristic Hill-type equation dependent on
remote insulin. This allows to define a unique time-varying glucose clearance for the three types of measured unlabeled and tracer glucose. In addition, the model provides parameter estimates for inhibition of glucose production by insulin, I50 and h, whose usefulness and reliability in providing information on glucose homeostasis still needs to be ascertained. Another difference with (Dalla Man et al., 2004) consists in the representation of glucose absorption rate, i.e. a smooth function with two parameters versus a piecewise linear function with several parameters and a fixed constraint on total glucose absorption during the test. The proposed simplified model imposes a constraint only on the shape of the absorption profile, e.g. it cannot represent biphasic gastro-intestinal glucose absorption. Nevertheless, thanks to the modeled between-subject variability of absorption parameters, α and β, the model is able to fit well also individual glucose profiles. Modelling glucose absorption with a parametric function with few degrees of freedom has the additional advantage of allowing estimation of all model parameters, at least with a population approach, without the need of fixing some of the model parameters. Parameter estimates obtained in healthy subjects are comparable with those reported in previous studies if normalisation is taken into account. For
instance, in (Dalla Man et al., 2004; Dalla Man et al., 2005) insulin sensitivity, SI , is normalised with respect to BW (average value 77 kg) and reported values range between 0.071 and 0.094 dl min−1 µU−1 ml. These values are within confidence limits of Table 1, as are the reported values between 112 and 123 dl for glucose distribution volume, VG . In contrast, remote insulin turnover, P2 , whose previous values range between 0.011 and 0.043 min−1 , match the present estimates only at their lower limit. Also, only the value of glucose effectiveness SG =0.0118 min−1 described in (Dalla Man et al., 2005) yields under standard assumptions GEZI values between 0.42 and 0.67 dl min−1 that are close to those of Table 1. Although not conclusive, these results thus support a cross-validation of different, yet comparable, experimental procedures and model-based data analysis approaches. However, the validity of the proposed model, especially as regards hepatic glucose production and effects of diabetes, cannot be fully established because of lack of similar studies performed in diabetic subjects. One confounding factor in comparing results from different studies depends also on different model parameterisations. In particular, applying the logtransformation to parameters in the NLME procedure corresponds to the implicit assumption of lognormal distribution for the original model parameters. This choice has two advantages: (i) parameters are automatically constrained to be nonnegative, and (ii) robustness with respect to outlying parameter values is improved, because the lognormal distribution can show heavier tails than the normal distribution, i.e. it can assign higher probabilities to large deviations from the mean for certain parameter values. Within a NLME population approach with log-normal parameter distribution, the presence of few outlying values can therefore be compensated mostly through a higher standard deviation of random effects rather than a bias on both standard deviation and mean, as usually seen under normality assumptions. Thus, even if one could compare population means of log-transformed with geometric means of un-transformed parameters, robustness with respect to outlier and efficiency of different identification procedures can strongly affect parameter estimates. Practical experience shows that population modelling with NLME can help solving model identification problems in small groups of subjects, which would be otherwise unfeasible at individual basis. The reason is that individual estimates benefit from group estimates on occasionally little informative experiments. Improved estimates of fixed and random effects could have been, obviously, obtained with a larger cohort. Experience teaches
that the NLME procedure is useful also with large number of subjects provided that the statistical model already adequately describes experimental data. In other words, once an appropriate model has been formulated using a limited number of subjects that characterise well the overall cohort, speed of convergence of NLME can even improve, with a reduced number of iterations, by increasing the size of the data set. The time-consuming step may then become the numerical integration of the model equations and their sensitivities with respect to parameters. However, an efficient implementation is possible through a straightforward parallelisation of the simulation of individual subjects and by exploiting modern multiprocessor architectures found nowadays even on laptop computers.
REFERENCES Abbate SL, Fujimoto WY, Brunzell JD, Kahn SE (1993). Effect of heparin on insulin-glucose interactions measured by the minimal model technique: implications for reproducibility using this method. Metabolism 42, 353–357. Bergman RN, Ider YZ, Bowden CR, Cobelli C (1979). Quantitative estimation of insulin sensitivity. Am J Physiol 236, E667–677. Dalla Man C, Caumo A, Cobelli C (2002). The oral glucose minimal model: estimation of insulin sensitivity from a meal test. IEEE Trans Biomed Eng 49, 419–429. Dalla Man C, Caumo A, Basu R, Rizza R, Toffolo G, Cobelli C (2004). Minimal model estimation of glucose absorption and insulin sensitivity from oral test: validation with a tracer method. Am J Physiol Endocrinol Metab 287, E637–643. Dalla Man C, Caumo A, Basu R, Rizza R, Toffolo G, Cobelli C (2005). Measurement of selective effect of insulin on glucose disposal from labeled glucose oral test minimal model. Am J Physiol Endocrinol Metab 289, E909-914. Davidian M, Giltinan DM (1995). Nonlinear Models for Repeated Measurement Data. Chapman & Hall. London. Hoenig M, Thomaseth K, Brandao J, Waldron M, Ferguson DC (in press). Assessment and mathematical modeling of glucose turnover and insulin sensitivity in lean and obese cats. Domest Anim Endocrinol. Pinheiro JC, Bates DM (2000). Mixed-Effects Models in S and S-plus. Springer. New York. Thomaseth K. Multidisciplinary modelling of biomedical systems. Comput Methods Programs Biomed 71:189-201, 2003.
CNR Institute of Biomedical Engineering, Padova, Italy ∗∗ CNR Institute of Clinical Physiology, Pisa, Italy ∗∗∗ Univ of Texas Health Science Center at San Antonio, USA
Abstract: A mathematical model based on minimal physiological assumptions for describing simultaneously multiple glucose measurements during modified oral glucose tolerance test is presented. Parameter identification was carried out using a population approach, which allowed precise characterisation of average glucose kinetic parameters in the studied cohort, as well as between–subject variability associated with glucose tolerance state and other covariates. The statistical and computational complexity added with the formulation of non-linear hierarchical population kinetic models is compensated by improved robustness, over the traditional two–stage approach, in selecting the most adequate model structure/order at the individual level and the most significant determinants of c 2006 IFAC glucose kinetics at the population level. Copyright Keywords: Mathematical models, Medical applications, Parameter identification, Physiological models, Reduced-order models, Statistical analysis
1. INTRODUCTION The oral glucose tolerance test (OGTT) is increasingly used as clinical test for studying glucose tolerance and insulin secretion in cohorts at risk of type 2 diabetes and for assessing the efficacy of new antidiabetic drugs. In its simplest form (oral load of 75 g glucose dissolved in water and a single measure of plasma glucose concentration after 2 hours) the OGTT is recommended by the World Health Organization for the diagnosis of diabetes in presence of elevated fasting plasma glucose. Given its moderate invasiveness for patients and usefulness in providing dynamic information on the kinetics of glucose and other endogenous substances, more complex OGTT-based protocols have been employed, e.g. with frequent blood sampling and administration of exogenous sub-
stances, combined with mathematical modelling of dynamic concentration profiles. In this study we consider the double-tracer labeled OGTT, which consists of a primedcontinuous infusion of tritiated glucose ([3-3 H]glucose) (to quantify endogenous glucose turnover by achieving steady state concentrations of [33 H]-glucose) and carbon labeled glucose ([1-14 C]glucose) added to the oral glucose load (to determine exogenous glucose absorption rates). The applicability of mathematical modelling to these kinds of experiments have been already demonstrated by others (Dalla Man et al., 2002; Dalla Man et al., 2004; Dalla Man et al., 2005). In particular, the so-called Bergman’s minimal model of glucose disappearance (MINMOD) (Bergman et al., 1979), which had been developed for modelling
frequently sampled IVGTT (FSIGT) data, has proven to adequately describe, after straightforward modifications, also glucose kinetics during OGTT. More specifically, by explicitly modelling the exogenous glucose appearance rate (Ra ) of orally administered unlabeled and [1-14 C]-glucose, it has been shown that MINMOD is able to accurately match the Ra profile determined with another, independent, modelling approach and to provide a figure of insulin sensitivity of glucose disappearance equivalent to that measured with the more laborious glucose clamp technique. Disadvantages of the model proposed in (Dalla Man et al., 2005) include the relatively high number of model parameters to be estimated from experimental data collected in a single individual and the use of two different models (similar structure but with different parameter values) for unlabeled and tracer glucose kinetics. Aims of this study were to re-assess the minimal modelling approach applied to double-tracer OGTT experiments within a unifying framework using: (i) a single model for describing simultaneously [3-3 H]-glucose, [1-14 C]-glucose and unlabeled glucose kinetics; (ii) a simplified model with reduced number of parameters, each having a unique value and physiological interpretation for different glucose measurements; (iii) a more robust population parameter estimation approach to better separate statistical variability of parameter estimates due to measurement and modelling errors from between-subject variability of kinetic parameters dependent from different metabolic states and other covariates, the characterisation of which is of major interest in population studies. The next section briefly summarises the studied cohort and the experimental protocol, and provides a detailed description of the mathematical model used to fit experimental data and the statistical population approach used to estimate model parameters. Section 3 presents results of the population modelling analysis, which are commented in section 4 together with a critical discussion about the proposed modelling approach.
2. MATERIALS AND METHODS 2.1 Subjects and protocol Twelve type 2 diabetic patients (7 males and 5 females; age = 53.6±2.5 years; body weight (BW) = 81.9±3.6 kg; body mass index (BMI) = 30.5±1.1 kg m−2 ) and ten normal glucose-tolerant subjects (5 males and 5 females; age = 39.6±3.7 years; BW = 90.1±3.9 kg; BMI = 31.1±0.9 kg m−2 ), matched for BMI and used as control subjects, underwent after overnight fast a doubletracer OGTT consisting of a primed-constant in-
fusion of [3-3 H]-glucose (0.25 µCi min−1 ) starting at time −120 min and an oral glucose load at time 0 of 75 g of glucose diluted in water containing 75 µCi [1-14 C]-glucose. Blood samples were collected every 15 minutes starting at time −30 min until 240 min for determination of plasma concentrations of unlabeled and labeled glucose and insulin. Concentration measurements are expressed as follows: glucose (mg dl−1 ), [3-3 H]glucose (dpm ml−1 ), [1-14 C]-glucose (dpm ml−1 ) and plasma insulin (µU ml−1 ).
2.2 Mathematical Modelling The model of glucose kinetics was formulated by assuming, for any of the three different glucose measurements, a one-compartment structure with fixed unknown glucose distribution volume (VG ). The generic mass balance equation is described by ˙ X(t) = −Rd (t) + Rp (t) + Ra (t)
(1)
where X(t) represents compartmental mass, expressed in (mg) and (dpm) for unlabeled and ˙ tracer glucose, and with time derivative X(t) (mg min−1 ) and (dpm min−1 ), respectively; glucose disappearance rate (Rd ) is defined as the product of instantaneous glucose clearance (CL) (dl min−1 ) and the corresponding plasma glucose concentration; the endogenous glucose production rate (Rp ) is nonzero only for unlabeled glucose; the exogenous glucose appearance rate (Ra ) is constant and known for [3-3 H]-glucose and timevarying and unknown for unlabeled and [1-14 C]glucose, yet with known total area under the curve (AUC) over an infinite horizon equal to the administered doses. Modelling glucose disappearance rate (Rd ). Given the above definition of Rd it is sufficient to define the mathematical structure for glucose clearance (CL). According to the minimal model of glucose disappearance (MINMOD) CL is described as the sum of a constant insulin-independent term, the glucose effectiveness at zero insulin (GEZI) (Abbate et al., 1993), and a time-varying term given as the product of insulin sensitivity index (SI ) (dl min−1 µU−1 ml) and insulin concentration in a remote compartment (IR (t)) CL(t) = GEZI + SI IR (t)
(2)
It is worth noting that the adopted parameterisation is different, but mathematically equivalent, to that conventionally used in the MINMOD literature, i.e. parameters expressing clearances (dl min−1 ) are obtained by multiplying the conventional counterpart of fractional clearances (min−1 ) by the glucose distribution volume (VG );
moreover the more commonly used glucose effectiveness at basal insulin (SG ) is related to the other model parameters through the equation SG = GEZI + SI Ib
(3)
The dynamics of insulin concentration in a remote compartment is described by I˙ R (t) = −P2 (IR (t) − I(t))
(4)
tion rate (Rp ) being dependent on deviations of remote insulin from basal through an empirical Hill-type equation. Given the glucose effectiveness at basal insulin (SG ), Eq. (3), and basal glucose concentration (Gb ), endogenous glucose production rate was modeled as Rp (t) = SG Gb
I50 h I50 h + (IR (t) − Ib )h
(8)
with initial condition equal to basal plasma insulin concentration (Ib ); plasma insulin concentration profile (I(t)) reconstructed by linear interpolation of measurements; and remote insulin turnover rate (P2 ). Thus IR (t) is simply a lagged profile of plasma insulin, and is related to the traditional MINMOD insulin action by: X(t) = SI (IR (t) − Ib ).
where I50 characterises the remote insulin concentration for 50% inhibition of basal glucose production, and the exponent h is an empirical parameter.
Modelling exogenous glucose appearance rate (Ra ). Exogenous unlabeled and [1-14 C]-glucose, administered through the oral load, share the same absorption time profile. A mathematical function, Ra (t, p), parametrised by p, is a suitable candidate for representing gastrointestinal glucose absorption if it is nonnegative and its total AUC is normalised to one, such that multiplication by administered dose yields effective absorption rate. For this purpose we adopted the following definition for Ra (t, p):
G ˙ G (t) = −CL(t) X (t) + RG X p (t) + Ra (t) (9) VG
T
T
β
Ra (t, p) dt = 1 − e−( α )
(5)
0
which ensures, if the two parameters α and β are positive, that the incremental AUC of glucose absorption strictly increases with time and tends exponentially to 1 for T → ∞. Given the oral unlabeled glucose dose (DG ) and the specific activity of administered [1-14 C]14 glucose (Sa C ), the exogenous glucose appearance rate for unlabeled and [1-14 C]-glucose are: G RG a (t) = D Ra (t, p) 14
C
14
Ra (t) = Sa
C
G
D Ra (t, p)
(6) (7)
Modelling endogenous glucose production rate (Rp ). In MINMOD endogenous glucose production is not modeled explicitly, but it is assumed that the control of net hepatic glucose balance (NHGB) parallels that of peripheral glucose utilisation. The MINMOD-derived insulin sensitivity index (SI ) for unlabeled glucose accounts thus also for the effects of insulin on NHGB. This leads to different values if SI is estimated from tracer glucose kinetics, which depends on glucose disposal alone. In view of our aim of defining a unique sensitivity index SI we assumed endogenous glucose produc-
Final model equations. The mathematical model equation used to fit double-tracer OGTT experiments can be summarised as follows G
14
˙ 14 C (t) = −CL(t) X X
14 (t) + Ra C (t) VG
3
˙ 3 H (t) = −CL(t) X X
C
H
(t)
VG
(10)
3
H + uinf
(11) 14
with initial conditions XG (0) = Gb VG , X C (0) = 3 3 0, X H (0) = GbH VG where CL(t) is derived from Eq. (2) with remote insulin determined through G Eq. (4); RG p (t) is expressed by Eq. (8); Ra (t) 14 3 C H and Ra (t) by Eq. (6) and (7), respectively; uinf represents the steady state constant infusion rate. The final set of estimated model parameters consist of: {P2 , GEZI, SI , VG , α, β, I50 , h}, with the first four parameters pertaining to the original MINMOD and two new parameters each for describing glucose absorption and production. Additional assigned model parameters are: test dose 14 3 H DG (75 g); individually measured Sa C and uinf ; 3 H basal concentration values Gb , Ib and Gb obtained for each individual as the average of measurements taken before the beginning of the test. Moreover, individual plasma insulin concentration profile (I(t)) is determined from insulin concentration measurements. Model equations were implemented using the modelling software Pansym (Thomaseth, 2003).
2.3 Parameter estimation Having defined the kinetic model for describing the concentration profile of unlabeled glucose, [33 H]-glucose and [1-14 C]-glucose, the i-th data set of the i-th subject can be formally represented through a non-linear regression model such as yi = f (ti , θi ) + ei
(12)
where yi is the vector of all measurements collected in the i-th subject, f (ti , θi ) represents the model predictions at the sampling times ti for given vector of parameters θi , and ei represents measurement (and modelling) noise, which was chosen among various alternatives as homoscedastic with a different variance estimated for each type of glucose measure. To guarantee the fulfillment of positivity constraints on parameters and to reduce the risk of identification problems in presence of large interindividual variability of kinetic parameters, logarithmic transformations were applied to model parameters before their estimation. Thus θ represent in the following the vector of log-transformed parameters, i.e. {log P2 , log GEZI, log SI , log VG , log α, log β, log I50 , log h}, which were the actually estimated parameters. Before their use in numerical simulation of model equations they were backtransformed through exponentiation. Results and statistical figures of precision (confidence intervals) will also be presented for the original, backtransformed, parameters. For obtaining information on the average kinetic response and on between-subject variability, indiˆ i , could have been vidual parameter estimates, θ obtained for each subject by non-linear weighted least squares fitting with subsequent statistical analysis to determine those covariates that affect kinetic parameters (this is the so-called twostep procedure). A more robust and dependable approach is the population modelling technique based on simultaneous analysis of all data using a statistical model that explicitly accounts for intra- and inter-individual variations of parameters (Davidian and Giltinan, 1995). In particular the non-linear mixed effects (NLME) approach was employed, which is based on maximum likelihood estimation of model parameters and covariance matrices of random variables under the assumption of normal distributions (Pinheiro and Bates, 2000). In this context, the k-th component of the individual parameter vector θi is in general described by a linear model as follows θik = θ0k + Xik ξ k + δik
(13)
where θ0k represents the average population value, Xik is a row vector of explicative variables (continuous and factorial covariates) to predict interindividual variability. Parameters θ0k and ξ k are called fixed effects because, by hypothesis, they don’t change between subjects of the studied population; δik are the so-called random effects that represent unpredictable between-subject variations. Yet, not all kinetic parameters were described using the complete structure of Eq. (13). To select the minimum number of estimated pa-
rameters able to quantify the statistically significant determinants of intra- and inter-individual variability, the inclusion of a covariate within Xik or the actual need for a random effect δik was assessed on the basis of physiological plausibility of results and statistical criteria for precision of parameter estimates and minimum model complexity (Davidian and Giltinan, 1995). For instance, the decision rule for including a covariate in the model was that the corresponding multiplicative regression parameter be estimated with a significance level of at least P=0.05 of being different from zero.
3. RESULTS The ability of the model to fit experimental data was excellent on pooled data (not shown) and very good on average data for normal control and diabetic groups, separately. Figure 1 shows the average concentrations (mean±SE) and model predictions for cold glucose, [3-3 H]-glucose, and [114 C]-glucose. For direct comparison each matching glucose measure is plotted with same vertical range but possibly with different offset. With regard to model performance at an individual level, the percentage error in model prediction (residual/fitted) had a low bias, being on average -0.02% for unlabeled glucose, 2.1% for [3-3 H]glucose, and 0.47% for [1-14 C]-glucose, and had an acceptable dispersion with standard deviations (a measure of mean coefficient of variation) of 6.2%, 16.4% and 7.8%, respectively. Fixed effects parameter estimates, characterizing both subgroups at the population level as one ensemble, are reported in Table 1 together with figures of precision, i.e. coefficient of variation and lower (L) and upper (U) limits of 95% statistical confidence intervals. It can be observed that parameter estimates have physiologically plausible values, such as GEZI = 0.72 dl min−1 that corresponds to a fractional glucose clearance of 0.0067 min−1 ; SI = 0.098 dl min−1 µU−1 ml in normal subjects, corresponds to a value of 9.17 10−4 min−1 µU−1 ml. The insulin sensitivity index in diabetic patients resulted in low values (only 4.9% of SI of normal subjects). In addition to the differences in SI between normal and diabetic subjects the only covariate found to significantly affect glucose kinetics was BW influencing GEZI. The value of %GEZI/kg reported in Table 1 quantifies the % variation in glucose effectiveness at zero insulin for increases in 1 kg BW. The reported reference value of GEZI corresponds to an arbitrarily chosen nominal BW of 85 kg. An interesting result is the quite low value estimated of 16 µU ml−1 for insulin concentration for 50% inhibition of basal glucose production (I50 ).
100
150
200
350 250
Glucose (mg/dl) 50
150
150
Glucose (mg/dl)
0 50 0
250
0
50
50
100
150
200
250
150
200
250
time(min)
200
250
200
250
6500
50
100
150
5000 2000 0
[1-14C]-Glucose (dpm/ml)
5000 2000
100
250
time(min)
0
[1-14C]-Glucose (dpm/ml)
50
200
5500 0
time(min)
0
150
4500
[3-3H]-Glucose (dpm/ml)
3500 0
100
time(min)
2500 1500
[3-3H]-Glucose (dpm/ml)
time(min)
0
50
100
150
time(min)
Fig. 1. Experimental data (mean±SE) and best average model predictions in normal controls (left panels) and diabetic subjects (right panels). Equivalent data are drawn with same scales but, whenever convenient, with different offsets. Table 1. Fixed effects population parameter estimates and precision. Parameter P2 GEZI %GEZI/kg SInormal SIdiabetic VG α β I50 h
Value 0.0163 0.72 2.09 0.0981 0.0048 107 123.1 1.451 16.11 1.425
CV% 28.9 9.74 0.70 14.8 18.2 18.1 9.3 3.1 24.7 22.7
95%(L) 0.0093 0.59 0.71 0.073 0.0033 75.1 102.5 1.36 9.92 0.91
95%(U) 0.0287 0.87 3.5 0.13 0.0068 152.5 147.8 1.55 26.2 2.22
4. DISCUSSION This is the first report, to our knowledge, on simultaneous modelling of multiple sampled glucose data from double-tracer oral glucose tolerance test in healthy and diabetic subjects using a population modelling approach. A similar approach has been adopted for the analysis of nonsteady state glucose kinetics during euglycemic hyperinsulinemic glucose clamp performed in cats (Hoenig et al., in press). The proposed model is closely related to previous work by others, e.g. (Dalla Man et al., 2004), the main difference being that glucose production is modeled here explicitly with a heuristic Hill-type equation dependent on
remote insulin. This allows to define a unique time-varying glucose clearance for the three types of measured unlabeled and tracer glucose. In addition, the model provides parameter estimates for inhibition of glucose production by insulin, I50 and h, whose usefulness and reliability in providing information on glucose homeostasis still needs to be ascertained. Another difference with (Dalla Man et al., 2004) consists in the representation of glucose absorption rate, i.e. a smooth function with two parameters versus a piecewise linear function with several parameters and a fixed constraint on total glucose absorption during the test. The proposed simplified model imposes a constraint only on the shape of the absorption profile, e.g. it cannot represent biphasic gastro-intestinal glucose absorption. Nevertheless, thanks to the modeled between-subject variability of absorption parameters, α and β, the model is able to fit well also individual glucose profiles. Modelling glucose absorption with a parametric function with few degrees of freedom has the additional advantage of allowing estimation of all model parameters, at least with a population approach, without the need of fixing some of the model parameters. Parameter estimates obtained in healthy subjects are comparable with those reported in previous studies if normalisation is taken into account. For
instance, in (Dalla Man et al., 2004; Dalla Man et al., 2005) insulin sensitivity, SI , is normalised with respect to BW (average value 77 kg) and reported values range between 0.071 and 0.094 dl min−1 µU−1 ml. These values are within confidence limits of Table 1, as are the reported values between 112 and 123 dl for glucose distribution volume, VG . In contrast, remote insulin turnover, P2 , whose previous values range between 0.011 and 0.043 min−1 , match the present estimates only at their lower limit. Also, only the value of glucose effectiveness SG =0.0118 min−1 described in (Dalla Man et al., 2005) yields under standard assumptions GEZI values between 0.42 and 0.67 dl min−1 that are close to those of Table 1. Although not conclusive, these results thus support a cross-validation of different, yet comparable, experimental procedures and model-based data analysis approaches. However, the validity of the proposed model, especially as regards hepatic glucose production and effects of diabetes, cannot be fully established because of lack of similar studies performed in diabetic subjects. One confounding factor in comparing results from different studies depends also on different model parameterisations. In particular, applying the logtransformation to parameters in the NLME procedure corresponds to the implicit assumption of lognormal distribution for the original model parameters. This choice has two advantages: (i) parameters are automatically constrained to be nonnegative, and (ii) robustness with respect to outlying parameter values is improved, because the lognormal distribution can show heavier tails than the normal distribution, i.e. it can assign higher probabilities to large deviations from the mean for certain parameter values. Within a NLME population approach with log-normal parameter distribution, the presence of few outlying values can therefore be compensated mostly through a higher standard deviation of random effects rather than a bias on both standard deviation and mean, as usually seen under normality assumptions. Thus, even if one could compare population means of log-transformed with geometric means of un-transformed parameters, robustness with respect to outlier and efficiency of different identification procedures can strongly affect parameter estimates. Practical experience shows that population modelling with NLME can help solving model identification problems in small groups of subjects, which would be otherwise unfeasible at individual basis. The reason is that individual estimates benefit from group estimates on occasionally little informative experiments. Improved estimates of fixed and random effects could have been, obviously, obtained with a larger cohort. Experience teaches
that the NLME procedure is useful also with large number of subjects provided that the statistical model already adequately describes experimental data. In other words, once an appropriate model has been formulated using a limited number of subjects that characterise well the overall cohort, speed of convergence of NLME can even improve, with a reduced number of iterations, by increasing the size of the data set. The time-consuming step may then become the numerical integration of the model equations and their sensitivities with respect to parameters. However, an efficient implementation is possible through a straightforward parallelisation of the simulation of individual subjects and by exploiting modern multiprocessor architectures found nowadays even on laptop computers.
REFERENCES Abbate SL, Fujimoto WY, Brunzell JD, Kahn SE (1993). Effect of heparin on insulin-glucose interactions measured by the minimal model technique: implications for reproducibility using this method. Metabolism 42, 353–357. Bergman RN, Ider YZ, Bowden CR, Cobelli C (1979). Quantitative estimation of insulin sensitivity. Am J Physiol 236, E667–677. Dalla Man C, Caumo A, Cobelli C (2002). The oral glucose minimal model: estimation of insulin sensitivity from a meal test. IEEE Trans Biomed Eng 49, 419–429. Dalla Man C, Caumo A, Basu R, Rizza R, Toffolo G, Cobelli C (2004). Minimal model estimation of glucose absorption and insulin sensitivity from oral test: validation with a tracer method. Am J Physiol Endocrinol Metab 287, E637–643. Dalla Man C, Caumo A, Basu R, Rizza R, Toffolo G, Cobelli C (2005). Measurement of selective effect of insulin on glucose disposal from labeled glucose oral test minimal model. Am J Physiol Endocrinol Metab 289, E909-914. Davidian M, Giltinan DM (1995). Nonlinear Models for Repeated Measurement Data. Chapman & Hall. London. Hoenig M, Thomaseth K, Brandao J, Waldron M, Ferguson DC (in press). Assessment and mathematical modeling of glucose turnover and insulin sensitivity in lean and obese cats. Domest Anim Endocrinol. Pinheiro JC, Bates DM (2000). Mixed-Effects Models in S and S-plus. Springer. New York. Thomaseth K. Multidisciplinary modelling of biomedical systems. Comput Methods Programs Biomed 71:189-201, 2003.