The iterative two-stage population approach to IVGTT minimal modeling: improved precision with reduced sampling

Paolo Vicini1 and Claudio Cobelli2

1 Department of Bioengineering, University of Washington, Seattle, Washington 98195; and 2 Department of Electronics and Informatics, University of Padova, Padua, Italy 35123


    ABSTRACT
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS
DISCUSSION
REFERENCES

The minimal model method is widely used to estimate glucose effectiveness (SG) and insulin sensitivity (SI) from intravenous glucose tolerance test (IVGTT) data. In the standard protocol (sIVGTT, 0.33 g/kg glucose bolus given at time 0), which allows the simultaneous assessment of beta -cell function, the precision of the individualized estimates often degrades and particularly so in the presence of reduced sampling schedules. Here, we investigated the use of a population approach, the iterative two-stage (ITS) approach, to analyze 16 sIVGTTs in healthy subjects and to obtain refined estimates of SG and SI in the population and in the individual subjects. The ITS is based on calculation of the population mean and standard deviation of the parameters at each iteration and then use of them as prior information for the individual analyses. Theoretically, the use of a prior in the ITS should improve the precision of the individual estimates. The customary approach (standard two stage, STS), where modeling is performed separately for each individual subject, does not take the population knowledge into account. We used both frequent (FSS, 30 samples) and (quasi-optimally) reduced (RSS, 14 samples) sampling schedules. For the FSS, STS gave estimates (mean ± SD) for SG = 2.66 ± 1.09 × 10-2 · min-1 and SI = 6.46 ± 6.99 10-4 · min-1 · µU-1 · ml, with an average precision of 51 (range 5-176) and 33% (3-91), respectively. RSS radically worsened the precision of both SG and SI. However, RSS and ITS gave SG = 2.59 ± 0.73 and SI = 6.06 ± 7.28, with an average precision of 23 (12-42) and 27% (), respectively. In conclusion, population minimal modeling of sIVGTT data improves the precision of individual estimates of glucose effectiveness and insulin sensitivity, as the theory predicts, and, even with reduced sampling, the improvement is substantial.

glucose effectiveness; insulin sensitivity; parameter estimation


    INTRODUCTION
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS
DISCUSSION
REFERENCES

THE MINIMAL MODEL of glucose disappearance (5) is widely used to estimate metabolic indexes of glucose effectiveness (SG) and insulin sensitivity (SI) in humans in both normal and physiopathological conditions. Since 1994, ~280 peer-reviewed scientific reports from the intermediary metabolism community have used the minimal model approach or have discussed related methodology, demonstrating substantial interest in the method; 65 of these reports have appeared in the past two years. Use of the minimal model for SI determination is often preferred over the more invasive and labor-intensive glucose clamp technique, especially in large clinical trials and epidemiological studies (26). Minimal model-based insulin sensitivity (SI) functions very well as a surrogate measure of cardiovascular disease risk in normal subjects (18) and in type 2 diabetics (17). Currently, several studies are using it to search for diabetes-relevant genes in the human genome (15, 16). Caumo et al. (7) have criticized this method because of its likely undermodeling of plasma glucose kinetics. However, its rising popularity among several circles warrants some effort to improve not only the data collection portion of the experiment but also the related procedures of data analysis, because physiologically significant conclusions will rest on both of these aspects.

Because of its simplicity, the standard intravenous glucose tolerance test (IVGTT) protocol (sIVGTT), based on the injection of a glucose bolus at time 0 and subsequent sampling for 3 or 4 h, is by far the most popular experimental protocol. However, the precision of SG and SI estimates is not always satisfactory, especially with reduced sampling schedules. The insulin-modified (mIVGTT) protocol (13, 37), which involves an infusion of insulin between 20 and 25 min after the glucose bolus, has allowed considerable improvement in the precision of the individual estimates of SG and SI. However, the modified protocol becomes substantially more complex than the standard procedure and may lead to hypoglycemic episodes in normal volunteers; also, it does not carry over easily to pediatric populations. In addition, the presence of exogenous insulin makes the determination of concomitant beta -cell activity difficult, and Pacini et al. (21) have recently suggested that, in studies comparing groups, use of the sIVGTT is more desirable. Improvement of the precision of the estimates in the sIVGTT is therefore paramount.

In the standard approach to minimal modeling, each subject is analyzed individually, and the additional information that all subjects belong to a homogeneous population is never exploited, even when available. With a single exception (10), minimal model analysis of large groups has so far neglected the fact that each individual belongs to a population of subjects who share certain quantitative traits. Nevertheless, the fact that all subjects in similarly performed studies have similar characteristics, like body weight, age, medications, etc., carries a certain amount of information. Thus the techniques of population analysis could bring a significant contribution to the minimal model method. Population analysis is the methodology used to quantify between-subject (also called intersubject in the literature) variability relative to a given population model. Its use is widespread in pharmacokinetic/pharmacodynamic studies, because it helps solve problems crucial in clinical studies such as sparse sampling and protocol deviation. Similar issues will become increasingly important in physiology as metabolic studies like the IVGTT move out of the investigative stage (small number of subjects) into the clinical arena (large number of subjects). The statistics literature has especially studied population analysis, developing the related methodologies of analysis of repeated measurement data, nonlinear mixed-effects modeling, and analysis of longitudinal data (9, 28).

In this work, we will describe the results of population analysis applied to glucose minimal modeling, and we will use the population results to provide Bayesian priors for the individual estimates of the metabolic indexes. The database consists of 16 sIVGTTs performed in normal subjects; we analyze these data with both a frequent (FSS) and a reduced sampling schedule (RSS). The use of an RSS is of interest in experimental design, because it potentially reduces experiment costs, laboratory worker biohazard, and patient discomfort. We will show that the use of a population analysis iterative algorithm, the iterative two-stage (ITS) approach, allows significant improvement in the precision of individual parameter estimates for each subject, even with RSS.


    MATERIALS AND METHODS
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS
DISCUSSION
REFERENCES

Data. The data consist of a series of sIVGTTs performed in a population of 16 young adults. Experiments were performed at the Washington University School of Medicine, St. Louis, MO, and the Department of Metabolic Diseases, University of Padova, Padua, Italy. We published these data previously in Vicini et al. (35), to which we refer the reader for experiment details. Briefly, a glucose bolus (between 300 and 330 mg/kg) was administered at time 0. Samples were drawn with an FSS of 30 samples, from which an optimal RSS of only 14 samples was extracted (0, 2, 3, 4, 5, 8, 10, 12, 14, 16, 18, 20, 24, 28, 32, 40, 45, 50, 60, 70, 80, 90, 100, 110, 120, 140, 160, 180, 210, 240 min; the RSS in italics). The RSS is based on optimal glucose sampling schedule studies (8, 25). Because different laboratories gathered the data, the sampling schedule differed slightly among subjects; the ones we report are only indicative.

Minimal model of glucose disappearance. We describe the minimal model of glucose disappearance (5) by
<A><AC>Q</AC><AC>˙</AC></A>(t)<IT>=</IT>−[S<SUB>G</SUB><IT>+</IT>X(t)]Q(t)<IT>+</IT>S<SUB>G</SUB>Q<SUB>b</SUB> G(<IT>0</IT>)<IT>=</IT>D<IT>/</IT>V<IT>+</IT>G<SUB>b</SUB>

<A><AC>X</AC><AC>˙</AC></A>(t)<IT>=</IT>−p<SUB><IT>2</IT></SUB>{X(t)<IT>−</IT>S<SUB>I</SUB>[I(t)<IT>−</IT>I<SUB>b</SUB>]} X(<IT>0</IT>)<IT>=0</IT> (1)

G(t)<IT>=</IT>Q(t)<IT>/</IT>V
where D is the glucose dose (mg/kg body wt), Q(t) (mg/kg) is glucose mass in plasma, G(t) (mg/dl) is glucose concentration, I(t) (µU/ml) is insulin concentration, Gb (= Qb/V) and Ib are their basal values, and X(t) is insulin action (min-1). Insulin concentration acts as a known input (without error) in the second equation, and we estimate model parameters by fitting the model response, G(t), to glucose concentration data. The model has four uniquely identifiable parameters: SG (min-1), SI (min-1 · µU-1 · ml), p2 (min-1), the insulin action parameter, and V (dl kg-1), the glucose distribution volume. SG and SI are the minimal model indexes of glucose effectiveness and insulin sensitivity, respectively, and reflect the effect of glucose and insulin on both glucose disposal and production. In particular, SG measures the ability of glucose per se, at basal insulin, to stimulate glucose disposal and to inhibit glucose production. Similarly, SI measures the ability of insulin to enhance the glucose per se stimulation of glucose disposal and the glucose per se inhibition of glucose production. The minimal model parameter vector is thus p = [SG, SI, p2, V].

Standard two-stage approach. The typical method of estimating the mean (first-order moment) and the variance (second-order moment) of the population distribution of SG and SI consists of estimating SG and SI separately in each subject via nonlinear regression and then calculating the sample mean (µ) and covariance (Sigma 2) of all the SG and SI estimates. The results are grouped according to the population (e.g., a population of normal controls and a population of diabetic patients should be analyzed separately). This method is often called the standard two-stage (STS) approach. The values of µ and Sigma 2 represent the mean and variance (or rather, the first- and second-order moments, if we cannot assume normality) of the population distribution. For the STS method, we used weighted nonlinear least squares, as implemented in the kinetic analysis software SAAM II (1), to determine each subject's parameters. We minimized the weighted residuals sum of squares with respect to the vector of model parameters for a subject j, pj
WRSS(p<SUB>j</SUB>)<IT>=</IT><LIM><OP>∑</OP><LL>i=<IT>1</IT></LL><UL>N<SUB>j</SUB></UL></LIM> <FR><NU>[G<SUP>OBS</SUP><SUB>i,j</SUB><IT>−</IT>G(<B>p</B><SUB>j</SUB><IT>, </IT>t<SUB>i,j</SUB>)]<SUP><IT>2</IT></SUP></NU><DE><IT>&sfgr;</IT><SUP><IT>2</IT></SUP><SUB>i,j</SUB></DE></FR> (2)
where Nj is the number of data points available for the jth subject, ti,j and Gi,jOBS are the ith time point and data point, respectively, of the jth subject, sigma 2i,j is the variance of the measurement error of the ith data point, and G(pj,ti,j) is the minimal model prediction of glucose concentration for a given pj.

We assumed measurement errors to be independent, Gaussian, zero mean, and with a standard deviation given by a constant coefficient of variation (CV) = 2% of the measured glucose concentration. We calculated the variance (precision), Vj, of the resulting estimate pj of pj from the inverse of the appropriate Fisher Information Matrix
V<SUP>STS</SUP><SUB>j</SUB><IT>=</IT><FENCE><FENCE><FR><NU><IT>∂</IT>G(<B>p</B><SUB>j</SUB><IT>, </IT>t<SUB>j</SUB>)</NU><DE><IT>∂</IT><B>p</B><SUB>j</SUB></DE></FR></FENCE><FENCE><AR><R><C><IT>&sfgr;</IT><SUP>−<IT>2</IT></SUP><SUB><IT>1,</IT>j</SUB></C><C><IT>0</IT></C><C><IT>0</IT></C></R><R><C><IT>0</IT></C><C><IT>⋱</IT></C><C><IT>0</IT></C></R><R><C><IT>0</IT></C><C><IT>0</IT></C><C><IT>&sfgr;</IT><SUP>−<IT>2</IT></SUP><SUB>N<SUB>j</SUB>,j</SUB></C></R></AR></FENCE><FENCE><FR><NU><IT>∂</IT>G(<B>p</B><SUB>j</SUB><IT>, </IT>t<SUB>j</SUB>)</NU><DE><IT>∂</IT><B>p</B><SUB>j</SUB></DE></FR></FENCE><SUP>T</SUP></FENCE><SUP><IT>−1</IT></SUP> (3)
where the superscript T indicates vector or matrix transpose, and [partial G/partial pj] is a 4 × N matrix that contains the derivatives of the model output at given times with respect to the parameters. We use this asymptotic approximation, which is known to give a lower bound of precision (as shown by the Cramèr-Rao theorem; see, e.g., Ref. 6, p. 200), both because it is widely believed to be a good index for "true" precision and for consistency with the usual applications of the minimal model method. To mitigate the error of the single compartment assumption, we did not use the early glucose data (from 0 to >= 6 min) for model identification. We calculated the basal concentrations of glucose and insulin (Gb and Ib) as the end-test basal concentrations, i.e., the last available measurement (usually at 240 min). Figure 1 shows the individual fits, all superimposed to the individual data, giving an idea of the degree of variability in the individual measurements. The variability is especially concentrated in the phase of acute response to the glucose bolus (between 20 and 70 min from the start of the experiment).


View larger version (23K):
[in this window]
[in a new window]
 
Fig. 1.   Plot of individual data points (black-lozenge ) and individual minimal model fits (---) for the frequent sampling schedule (FSS) and standard two-stage (STS) analysis.

We calculated the population mean for each parameter as the sample mean of all the individual parameter estimates (pj)
&mgr;<SUB>STS</SUB><IT>=</IT><FR><NU><IT>1</IT></NU><DE>N</DE></FR> <LIM><OP>∑</OP><LL>j=<IT>1</IT></LL><UL>N</UL></LIM> <B><A><AC>p</AC><AC>ˆ</AC></A></B><SUB>j</SUB> (4)
where N is the number of subjects, and we calculated the population variance as the corresponding sample variance
&Sgr;<SUP><IT>2</IT></SUP><SUB>STS</SUB><IT>=</IT><FR><NU><IT>1</IT></NU><DE>N<IT>−1</IT></DE></FR> <LIM><OP>∑</OP><LL>j=<IT>1</IT></LL><UL>N</UL></LIM> (<B><A><AC>p</AC><AC>ˆ</AC></A></B><SUB>j</SUB><IT>−&mgr;</IT><SUB>STS</SUB>)(<B><A><AC>p</AC><AC>ˆ</AC></A></B><SUB>j</SUB><IT>−&mgr;</IT><SUB>STS</SUB>)<SUP>T</SUP> (5)
From a methodological standpoint, this approach is undesirable for a number of reasons. It neglects the fact that the precision of the estimates of SG and SI can differ substantially between individuals, and it therefore pools "good" individual estimates together with "poor" individual estimates in Eqs. 4 and 5. The approach is also well known to overestimate, to a varying degree dependent on the magnitudes of Vj, the true population variance (9).

ITS approach. The ITS is a parametric, iterative population analysis method based on the concepts of population prior knowledge and maximum a posteriori (MAP) probability empirical Bayes estimation. Steimer et al. (34) proposed it as a computationally attractive alternative to the nonlinear mixed-effects model approach (9). The steps of the ITS follow.

Step 1: initialization. Estimate the STS mean and covariance as in Eqs. 3 and 4. Define µ(0) = µSTS, Sigma 2(0) = Sigma STS2.

Step 2: k + 1, k >=  1. Perform parameter estimation on each subject j again, this time minimizing the following extended MAP Bayesian objective function with respect to pj
MAP(p<SUB>j</SUB>)<IT>=</IT><LIM><OP>∑</OP><LL>i=<IT>1</IT></LL><UL>N<SUB>j</SUB></UL></LIM> <FR><NU>[G<SUP>OBS</SUP><SUB>i,j</SUB><IT>−</IT>G(<B>p</B><SUB>j</SUB><IT>, </IT>t<SUB>i,j</SUB>)]<SUP><IT>2</IT></SUP></NU><DE><IT>&sfgr;</IT><SUP><IT>2</IT></SUP><SUB>i,j</SUB></DE></FR><IT>+</IT><LIM><OP>∑</OP><LL>i=<IT>1</IT></LL><UL>N<SUB>p</SUB></UL></LIM> <FR><NU>[<IT>&mgr;</IT><SUB>i</SUB>(k)<IT>−</IT><B>p</B><SUB>j,i</SUB>]<SUP><IT>2</IT></SUP></NU><DE><IT>&Sgr;</IT><SUP><IT>2</IT></SUP><SUB>i,i</SUB>(k)</DE></FR> (6)
where the distance of the current parameter estimate from the population mean (the prior) is also penalized; we denote with pj,i the ith element of the parameter vector for subject j and with Np the number of elements in the vector pj. The estimate pj obtained by minimizing this objective function is often called post hoc, or empirical Bayes, estimate. We can again calculate Vj, the precision of the parameter estimate pj, from the Fisher Information Matrix (see below); µi(k) is the value of the population mean at the kth iteration of the method, and Sigma i,i(k) is the ith diagonal element of the population covariance matrix at the kth iteration. Calculate, then, the updated population mean of the parameter vector
&mgr;(k<IT>+1</IT>)<IT>=</IT><FR><NU><IT>1</IT></NU><DE>N</DE></FR> <LIM><OP>∑</OP><LL>j=<IT>1</IT></LL><UL>N</UL></LIM> <B><A><AC>p</AC><AC>ˆ</AC></A></B><SUB>j</SUB>(k) (7)
and the covariance
&Sgr;<SUP>2</SUP>(k<IT>+1</IT>)<IT>=</IT><FR><NU><IT>1</IT></NU><DE>N</DE></FR> <LIM><OP>∑</OP><LL>j=<IT>1</IT></LL><UL>N</UL></LIM> {V<SUB>j</SUB>(k) (8)

<IT>+</IT>[<B><A><AC>p</AC><AC>ˆ</AC></A></B><SUB>j</SUB>(k)<IT>−&mgr;</IT>(k<IT>+1</IT>)][<B><A><AC>p</AC><AC>ˆ</AC></A></B><SUB>j</SUB>(k)<IT>−&mgr;</IT>(k<IT>+1</IT>)]<SUP>T</SUP>}
Step 3: k + 2, k >=  1. Check for convergence of the population mean, the population variance, and the individual parameter estimates; namely, determine whether or not the current and the previous estimate differ by <1%. If so, stop; if not, return to step 2.

The iterative nature of the algorithm is apparent; it is also apparent that the objective function in Eq. 5 takes explicitly into account the available population information and that the computation of the population covariance matrix in Eq. 8 includes available information on the precision of each individual estimate. The STS, on the other hand, ignores subjects other than the current one; it is essentially the initialization step of the ITS. The ITS has had limited use in the pharmacokinetic literature to estimate population parameter means and variances from reduced data sets (11). The ITS has been implemented in a prototype application built on SAAM II, with the maximum number of iterations set to 50. We observed essentially no change in the individual and population parameter values and precisions after 10-45 iterations of the method.

We computed the precisions of the individual parameter estimates from the inverse of the Fisher Information Matrix again. Unlike STS, which uses only the kinetic glucose measurements, the ITS covariance formula also uses the population information (20)
V<SUP>ITS</SUP><SUB>j</SUB><IT>=</IT><FENCE><FENCE><FR><NU><IT>∂</IT>G(<B>p</B><SUB>j</SUB><IT>, </IT>t<SUB>j</SUB>)</NU><DE><IT>∂</IT><B>p</B><SUB>j</SUB></DE></FR></FENCE><FENCE><AR><R><C><IT>&sfgr;</IT><SUP>−<IT>2</IT></SUP><SUB><IT>1,</IT>j</SUB></C><C><IT>0</IT></C><C><IT>0</IT></C></R><R><C><IT>0</IT></C><C><IT>⋱</IT></C><C><IT>0</IT></C></R><R><C><IT>0</IT></C><C><IT>0</IT></C><C><IT>&sfgr;</IT><SUP>−<IT>2</IT></SUP><SUB>N<SUB>j</SUB>,j</SUB></C></R></AR></FENCE></FENCE> (9)

<FENCE><IT>× </IT><FENCE><FR><NU><IT>∂</IT>G(<B>p</B><SUB>j</SUB><IT>, </IT>t<SUB>j</SUB>)</NU><DE><IT>∂</IT><B>p</B><SUB>j</SUB></DE></FR></FENCE><SUP>T</SUP><IT>+&Sgr;</IT><SUP>−<IT>2</IT></SUP></FENCE><SUP>−<IT>1</IT></SUP>

where Sigma 2 is the final estimate of the population covariance from the ITS. Although it is apparent that Eq. 9 will always improve precision compared with Eq. 3 (because Sigma 2 is a positive definite matrix by construction), the question arises of how to choose a proper value for Sigma 2. The ITS provides an answer to that problem.

Statistical methods. Because of both the nonnormality of the sample and the heterogeneity of variances between the FSS and RSS groups (see Ref. 30, p. 296), we compared the SI and SG values obtained from different estimators and different sampling schedules by means of the Wilcoxon signed-ranks test for matched pairs (Ref. 30, p. 128). We set the significance level at P = 0.05, and we used the nonparametric Spearman rank correlation to measure agreement.


    RESULTS
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS
DISCUSSION
REFERENCES

We report here the STS and the ITS results for both sampling schedules (FSS and RSS). Table 1 (FSS by use of STS and ITS) and Table 2 (RSS by use of STS and ITS) show individual results for all of the subjects for all of the minimal model parameters, whereas Fig. 2 displays a visual summary of our findings on the percent precisions of parameter estimates.

                              
View this table:
[in this window]
[in a new window]
 
Table 1.   STS and ITS results for FSS sIVGTT in normal subjects


                              
View this table:
[in this window]
[in a new window]
 
Table 2.   STS and ITS results for RSS sIVGTT in normal subjects



View larger version (20K):
[in this window]
[in a new window]
 
Fig. 2.   Precision of glucose effectiveness (SG) and insulin sensitivity (SI) estimates with different sampling schedules [FSS or reduced sampling schedule (RSS)] and population analysis approaches [STS and iterative two stage (ITS)]. Mean values are shown.

We first applied the STS method to both FSS and RSS data. We detected a significant difference (P < 0.05) between SI determined with STS analysis of the FSS compared with the RSS data. To investigate the effect of the number of samples on population analysis, we used linear regression to compare the individual estimates across sampling schedules. From Fig. 3 we can infer that SI is fairly stable (except when its value is very low, as in subject 11), whereas SG is slightly less stable. With the FSS data, the average precision of the estimates, calculated as the mean of all individual precisions, was 51, 33, 73, and 7%, respectively. However, the average precision of the estimates worsened considerably with the RSS data (79, 91, 227, and 11%, respectively). Therefore, the RSS considerably hampers individual parameter reliability.


View larger version (15K):
[in this window]
[in a new window]
 
Fig. 3.   Impact of the RSS on the SG and SI estimates from the standard intravenous glucose tolerance test (sIVGTT). Estimates are compared with the corresponding estimates obtained under FSS. R2 values reported are Spearman rank correlation results.

We then analyzed FSS and RSS data with the ITS method. Results for FSS data were comparable to the STS results. The population spread was slightly smaller than the STS estimate, as expected from theory, i.e., the STS overestimates the population covariance. Interestingly, we found no statistically significant difference between STS-FSS and ITS-RSS insulin sensitivity estimates, showing that the RSS is not robust unless the data are analyzed with a population analysis method. We can speculate that, in a data-rich situation, STS and ITS would perform similarly; however, the average precision of the estimates improved dramatically with ITS (24, 18, 35, and 4%, respectively). This approach maximizes the information available from both the rich data set of the FSS and the knowledge of the population mean. The ITS, however, reveals its power in a (relatively) data-poor situation like the RSS. The average precision of the estimates in this case was very good: 23, 27, 37, and 4%, respectively. These numbers are similar to those obtained with ITS-FSS, but the process employed fewer data points. The gain of the ITS with respect to the STS, with both sampling schedules, is quite evident from the summary of the estimated parameter precision in all four cases shown in Fig. 2.

Last, we investigated the performance of the estimation methods applied to the same sampling schedule. To give an idea of how individual estimates from the STS map into individual estimates from the ITS, we compared the individual STS and ITS parameter estimates via linear regression; rank correlation of SG and SI was always between 0.97 and 1.00.


    DISCUSSION
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS
DISCUSSION
REFERENCES

Population analysis methods find their natural application in the analysis of data-poor clinical studies, as when the number of samples available for each individual is rather small (e.g., 3 or 4). The natural arena is that of pharmacokinetics and pharmacodynamics, but there have been some recent contributions in the metabolism arena (24). We have found here that both rough (like the STS) and somewhat more sophisticated (like the ITS) methods give similar results for the population mean and variance. However, the population-based approach obtained the same results for population mean and variance with one-half of the total number of blood samples (an attractive feature). Moreover, population analysis allows a significant gain in the individual precisions of SG and SI (Fig. 2), thus giving more confidence in the overall results. A population method like the ITS allows the calculation of much more precise individual estimates for all of the subjects, even for those where the STS parameter precisions become unacceptable. For instance, our results demonstrate that an RSS is feasible for the sIVGTT, provided that a population-based method is used for data analysis; in fact, we found a statistically significant difference between the FSS and RSS SI estimates with STS analysis. Note that the overestimation error associated with the STS estimate of between-individual variability increases as within-individual variability (related to the precision of parameter estimates and also referred to as intraindividual variability) increases (which is particularly evident when comparing the estimate of Sigma  for p2 in Table 2, STS vs. ITS). This is the reason for the failure of STS when applied to the RSS case. Therefore, in an experimental protocol with very precisely estimated individual parameters, as in tracer studies (8), or with the mIVGTT (36), the STS would give quite reliable estimates of the population mean and variance, because within-individual variability would have only a small confounding effect. In contrast, in a situation where within-individual variability is as large as or larger than between-individual variability (i.e., the parameters are poorly estimated), the two effects would confound each other. Thus the sIVGTT protocol with RSS benefits considerably from population analysis.

The minimal model community widely uses measured values instead of predicted values when calculating sigma 2i,j in Eq. 2. This practice, however, introduces measurement error in the weights. Inspection of weighted residuals from the individual fits indicates that model misspecification is small (data not shown), thus suggesting that the difference would not be significant. Nevertheless, an extended least squares (ELS) or maximum likelihood criterion, where a term in the objective function with the logarithm of the data variance balances the use of predictions in the weighting, would better account for the variance in the data (2, 3)
ELS(p<SUB>j</SUB>)<IT>=</IT><LIM><OP>∑</OP><LL>i=<IT>1</IT></LL><UL>N<SUB>j</SUB></UL></LIM> log<IT> &sfgr;<SUP>2</SUP></IT>(<B>p</B><SUB>j</SUB><IT>, </IT>t<SUB>i,j</SUB>)<IT>+</IT><LIM><OP>∑</OP><LL>i=<IT>1</IT></LL><UL>N<SUB>j</SUB></UL></LIM> <FR><NU>[G<SUP>OBS</SUP><SUB>i,j</SUB><IT>−</IT>G(<B>p</B><SUB>j</SUB><IT>, </IT>t<SUB>i,j</SUB>)]<SUP><IT>2</IT></SUP></NU><DE><IT>&sfgr;<SUP>2</SUP></IT>(<B>p</B><SUB>j</SUB><IT>, </IT>t<SUB>i,j</SUB>)</DE></FR> (10)
Equation 10 gives twice the negative logarithm of the likelihood function of the data. If sigma 2(pj,tij) did not depend on the parameters and were constant across time (constant standard deviation), then the two methods (weighted and extended least squares) would be entirely equivalent.

Theory indicates that estimation of parameters with empirical Bayes estimation will result in greater precision (as calculated from the Fisher Information Matrix in Eqs. 3 and 9) than with the STS approach. However, the purpose of this study was to ascertain the magnitude and possible significance of such an improvement a priori unknown. Indeed, with extremely large variabilities of the metabolic indexes, we could expect negligible improvement. Also, we should note the (historical) fact that, in the clinic and even in large epidemiological studies (19), the STS is invariably used for SG and SI estimation, with both reduced and intensive sampling schedule (18-30 blood samples). We hope that this work can begin to introduce this large community to the concept that population-based methods and/or empirical Bayes approaches improve individual/group parameter quantification in the presence of large epidemiological data, and that it is feasible to use priors (when available and appropriate) to estimate individual patients' parameters (a practice sometimes referred to in pharmacokinetics as "Bayesian forecasting").

The improved precision of parameter estimates in this subject sample is due, at least in part, to the "shrinkage" of the individual estimates around the population mean (see Eq. 6). The not-so-hidden assumption is that a population mean indeed exists and that it is appropriate for all of the subjects involved in the study. Several levels of assumptions are embedded in kinetic modeling of population data: 1) that the mathematical model of the concentration time course is true; 2) that a certain measurement-error statistical structure exists (which allows the application of meaningful nonlinear regression); and 3) that, in the case of our version of the ITS, the shape of the population distribution of the parameters is multivariate normal. Examination of the consequences of all of these assumptions does not concern us here, since the present report is largely a "proof of concept" of the usefulness of population analysis in a reduced-data situation, and the ITS is only one possible algorithm for this purpose. We are thus not concerned here with either model accuracy or methods other than the ITS. However, some of the aforementioned assumptions, such as normality or unimodality of the population distribution of the parameters, can be relaxed, e.g., via mixture modeling, explicitly including covariants in the model, or performing nonparametric population analysis (see, e.g., Chapter 7 in Ref. 9). Davidian and Giltinan (9) comprehensively reviewed several different proposed estimation methods for population analysis of nonlinear models, all based on different degrees of approximation to the maximum-likelihood objective function. The attractiveness of the ITS derives from the need for a simple computational machinery based on successive iterations of the single-subject identification algorithm via a recursion formula. The ITS method is essentially an expectation-maximization (EM) algorithm (27). A related EM-type algorithm is the global two stage (GTS), which applies an iterative algorithm directly to estimated parameters and their within-individual covariances. Recently, Patron-Bizet et al. (22) compared the STS and the GTS used with a pharmacodynamic model and found that the GTS was superior. Steimer et al. (34) and Racine-Poon and Smith (23) compared the GTS and ITS, finding that they perform quite similarly. Although the ITS is computationally more intensive than the GTS, we chose to use the ITS here, because it allows a natural introduction of the concept of empirical Bayes individual estimation, which is at the root of the improved (asymptotic) precision. We have found the ITS method quite robust with respect to different starting values, with the main effect being the speed of convergence. After a few iterations, the ITS algorithm gets reasonably near its limit point. The literature suggests checking for convergence by calculation of the maximum-likelihood objective function, or an approximation thereof, at the point estimates (14).

Other available algorithms for the calculation of population parameters belong to the family of mixed-effects modeling, where the mean and the covariance of the population are fitted parameters in a maximum-likelihood context, and the computationally demanding maximum-likelihood integral is calculated via various approximations to the model function (29). Indeed, given that we have chosen a homogeneous population, mixed-effects modeling might be the best choice. Some of these approximations are currently available in the software NONMEM (28). Indeed, these algorithms are often preferred to the ITS, because they can be used even when individual estimates are not available in some subjects. They are based on the minimization of a well-defined objective function and return the precision of the estimates of population mean and covariance (and, in principle, also of each individual estimate). In contrast, the ITS algorithm does not readily allow calculation of the precision of population mean and variance. Because we were primarily interested in gauging the quantitative improvement of the individual posterior estimates, the ITS adequately met our purposes as a population analysis approach. Interestingly, De Gaetano et al. (10) used the first-order approximation to the population problem implemented in NONMEM to estimate minimal model population parameters in a group of 20 normal subjects by use of sIVGTT data. In that article, the precision of the population mean and covariance estimates improved with respect to the STS method. The mean and covariance estimates, however, did not change between NONMEM and STS (again, because the frequently sampled sIVGTT is a data-rich experiment). However, the authors did not report the values and precision of the individual parameter estimates (possibly because they are not reported by the NONMEM software).

Allowing for precise estimates of the parameters at the individual level in a reduced-data situation, and thus improving the a posteriori identifiability of the model parameters in a population, is only one of the possible applications of empirical Bayes estimation and population modeling in glucose metabolism and in intermediary metabolism in general. For example, there is ample pharmacokinetic/pharmacodynamic literature about the incorporation of demographic covariants into the mathematical model (12). Thus it is conceivable that population modeling of metabolic studies would, for instance, permit understanding the relationships of important features like visceral adiposity, weight, and age with metabolic indexes of insulin sensitivity and glucose effectiveness.


    ACKNOWLEDGEMENTS

We wish to thank our colleagues, Alan Schumitzky and Bradley M. Bell, for useful discussions and to gratefully acknowledge the editorial help of Eileen Thorsos and the useful suggestions of the anonymous referees.


    FOOTNOTES

This work was partially supported by National Institutes of Health Grants NCRR-12609 ("Resource Facility for Population Kinetics") and GM-53930. Preliminary results from the material in this work were presented at the Mathematical Modeling in Experimental Nutrition Conference held at the University of California, Davis, CA, on August 17-20, 1997, and at the American Diabetes Association 59th Scientific Sessions held in San Diego, CA, on June 18-22, 1999.

Address for reprint requests and other correspondence: P. Vicini, Dept. of Bioengineering, Box 352255, Univ. of Washington, Seattle, WA 98195-2255 (E-mail: vicini{at}u.washington.edu).

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked "advertisement" in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

Received 23 August 1999; accepted in final form 13 September 2000.


    REFERENCES
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS
DISCUSSION
REFERENCES

1.   Barrett, PHR, Bell BM, Cobelli C, Golde H, Schumitzky A, Vicini P, and Foster DM. SAAM II: simulation, analysis, and modeling software for tracer and pharmacokinetic studies. Metabolism 47: 484-492, 1998[ISI][Medline].

2.   Beal, SL, and Sheiner LB. Estimating population kinetics. Crit Rev Biomed Eng 8: 195-222, 1982[Medline].

3.   Beal, SL, and Sheiner LB. NONMEM User Guide. San Francisco, CA: NONMEM Project Group, Univ. of California San Francisco, 1992.

4.   Bergman, RN, Finegood DT, and Ader M. Assessment of insulin sensitivity in vivo. Endocr Rev 6: 45-86, 1985[ISI][Medline].

5.   Bergman, RN, Ider YZ, Bowden CR, and Cobelli C. Quantitative estimation of insulin sensitivity. Am J Physiol Endocrinol Metab Gastrointest Physiol 236: E667-E677, 1979[Abstract/Free Full Text].

6.   Carson, ER, Cobelli C, and Finkelstein L. The Mathematical Modeling of Metabolic and Endocrine Systems. New York: Wiley, 1983.

7.   Caumo, A, Vicini P, Zachwieja JJ, Avogaro A, Yarasheski K, Bier DM, and Cobelli C. Undermodeling affects minimal model indexes: insights from a two-compartment model. Am J Physiol Endocrinol Metab 276: E1171-E1193, 1999[Abstract/Free Full Text].

8.   Cobelli, C, and Ruggeri A. A reduced sampling schedule for estimating the parameters of the glucose minimal model from a labeled IVGTT. IEEE Trans Biomed Eng 38: 1023-1029, 1991[ISI][Medline].

9.   Davidian, M, and Giltinan D. Nonlinear Models for Repeated Measurement Data. New York: Chapman and Hall, 1995.

10.   De Gaetano, A, Mingrone G, and Castagneto M. NONMEM improves group parameter estimation for the minimal model of glucose kinetics. Am J Physiol Endocrinol Metab 34: E932-E937, 1996.

11.   Drusano, GL, Forrest A, Snyder MJ, and Reed MD. An evaluation of optimal sampling strategy and adaptive study design. Clin Pharmacol Ther 44: 232-238, 1988[ISI][Medline].

12.   Ette, EI, and Ludden TM. Population pharmacokinetic modeling: the importance of informative graphics. Pharm Res 12: 1845-1855, 1995[ISI][Medline].

13.   Finegood, DT, Hramiak IM, and Duprè J. A modified protocol for estimation of insulin sensitivity with the minimal model of glucose kinetics in patients with insulin dependent diabetes. J Clin Endocrinol Metab 70: 1538-1548, 1990[Abstract].

14.   Forrest, A, Ballow CH, Nix DE, Birmingham MC, and Schentag JJ. Development of a population pharmacokinetic model and optimal sampling strategies for intravenous ciprofloxacin. Antimicrob Agents Chemother 37: 1065-1072, 1993[Abstract].

15.   Ghosh, S, Hauser ER, Magnuson VL, Valle T, Ally DS, Karanjawala ZE, Rayman JB, Knapp JI, Musick A, Tannenbaum J, Te C, Eldridge W, Shapiro S, Musick T, Martin C, So A, Witt A, Harvan JB, Watanabe RM, Hagopian W, Eriksson J, Nylund SJ, Kohtamaki K, Tuomilehto-Wolf E, Toivanen L, Vidgren G, Enholm C, Bergman RN, Tuohmilehto J, Collins FS, and Boehnke M. A large sample of Finnish diabetic sib-pairs reveals no evidence for a non-insulin-dependent diabetes mellitus susceptibility locus at 2qter. J Clin Invest 102: 704-709, 1998[Abstract/Free Full Text].

16.   Ghosh, S, Langefeld CD, Ally D, Watanabe RM, Hauser ER, Magnuson VL, Nylund SJ, Valle T, Eriksson J, Bergman RN, Tuomilehto J, Collins FS, and Boehnke M. The W64R variant of the beta 3-adrenergic receptor is not associated with type II diabetes or obesity in a large Finnish sample. Diabetologia 42: 238-244, 1999[ISI][Medline].

17.   Haffner, SM, D'Agostino RB, Jr, Mykkanen L, Tracy R, Howard B, Rewers M, Selby J, Savage PJ, and Saad MF. Insulin sensitivity in subjects with type 2 diabetes: relationship to cardiovascular risk factors: the Insulin Resistance Atherosclerosis Study. Diabetes Care 22: 562-568, 1999[Abstract].

18.   Howard, G, Bergman RN, Wagenknecht LE, Haffner SM, Savage PJ, Saad MF, Laws A, and D'Agostino RB, Jr. Ability of alternative indices of insulin sensitivity to predict cardiovascular risk: comparison with the "minimal model". Insulin Resistance Atherosclerosis Study (IRAS) Investigators. Ann Epidemiol 8: 358-369, 1998[ISI][Medline].

19.   Karter, AJ, Mayer-Davis EJ, Selby JV, D'Agostino RB, Jr, Haffner SM, Sholinski P, Bergman R, Saad MF, and Hamman RF. Insulin sensitivity and abdominal obesity in African-American, Hispanic, and non-Hispanic white men and women. The Insulin Resistance and Atherosclerosis Study. Diabetes 45: 1547-1555, 1996[Abstract].

20.   Merlè, Y, and Mentrè F. Bayesian design criteria: computation, comparison, and application to a pharmacokinetic and a pharmacodynamic model. J Pharmacokinet Biopharm 23: 101-125, 1995[ISI][Medline].

21.   Pacini, G, Tonolo G, Sambataro M, Maioli M, Ciccarese M, Brocco E, Avogaro A, and Nosadini R. Insulin sensitivity and glucose effectiveness: minimal model analysis of standard and insulin-modified FSIGT. Am J Physiol Endocrinol Metab 274: E592-E599, 1998[Abstract/Free Full Text].

22.   Patron-Bizet, F, Mentre F, Genton M, Thomas-Haimez C, and Maccario J. Assessment of the global two-stage method to EC50 determination. J Pharmacol Toxicol Methods 39: 103-108, 1998[ISI][Medline].

23.   Racine-Poon, A, and Smith AFM A Bayesian approach to nonlinear random effects models. Biometrics 41: 1015-1023, 1985[ISI][Medline].

24.   Rostami-Hodjegan, A, Peacey SR, George E, Heller SR, and Tucker GT. Population-based modeling to demonstrate extrapancreatic effects of tolbutamide. Am J Physiol Endocrinol Metab 274: E758-E771, 1998[Abstract/Free Full Text].

25.   Ruggeri, A, and Cobelli C. Optimal sampling schedule design may reveal inadequacy of model structure: a case study on the minimal model of glucose disappearance. In: Modelling and Control in Biomedical Systems, edited by Cobelli C, and Mariani L. Oxford, UK: Pergamon, 1988, p. 121-125.

26.   Saad, MF, Anderson RL, Laws A, Watanabe RM, Kades WW, Chen Ida Y-D, Sands RE, Pei D, Savage PJ, and and Bergman RN for the IRAS A comparison between the minimal model and the glucose clamp in the assessment of insulin sensitivity across the spectrum of glucose tolerance. Diabetes 43: 1114-1121, 1994[Abstract].

27.   Schumitzky, A. EM algorithms and two stage methods in pharmacokinetic population analysis. In: Advanced Methods of Pharmacokinetic and Pharmacodynamic Systems Analysis, edited by D'Argenio DZ. New York: Plenum, 1995, p. 145-160.

28.   Sheiner, LB, and Beal SL. Evaluation of methods for estimating population pharmacokinetic parameters. II. Biexponential model and experimental pharmacokinetic data. J Pharmacokinet Biopharm 9: 635-651, 1981[ISI][Medline].

29.   Sheiner, LB, Rosenberg B, and Marathe VV. Estimation of population characteristics of pharmacokinetic parameters from routine clinical data. J Pharmacokinet Biopharm 5: 445-479, 1977[ISI][Medline].

30.   Snedecor, GW, and Cochran WG. Statistical Methods (6th ed.). Ames: Iowa State Univ, 1974.

32.   Steil, GM. Choosing an optimal sampling schedule for the identification of minimal model insulin sensitivity. In: The Minimal Model Approach and Determinants of Glucose Tolerance, edited by Bergman RN, and Lovejoy JC. Baton Rouge: Louisiana State Univ, 1997, p. 123-141.

33.   Steil, GM, Volund A, Kahn SE, and Bergman RN. Reduced sample number for calculation of insulin sensitivity and glucose effectiveness from the minimal model: suitability for use in population studies. Diabetes 42: 250-256, 1992[Abstract].

34.   Steimer, J-L, Mallet A, Golmard J-L, and Boisvieux J-F. Alternative approaches to estimation of population pharmacokinetic parameters: comparison with the nonlinear mixed-effect model. Drug Metab Rev 15: 265-292, 1984[ISI][Medline].

35.   Vicini, P, Caumo A, and Cobelli C. The hot IVGTT two-compartment minimal model: indexes of glucose effectiveness and insulin sensitivity. Am J Physiol Endocrinol Metab 273: E1024-E1031, 1997[ISI][Medline].

36.   Vicini, P, and Cobelli C. Population minimal modeling improves precision of the standard IVGTT estimates of insulin sensitivity and glucose effectiveness (Abstract). Diabetes 48, Suppl 1: 1297, 1999.

37.   Yang, YJ, Youn JH, and Bergman RN. Modified protocols improve insulin sensitivity estimation using the minimal model. Am J Physiol Endocrinol Metab 253: E595-E602, 1987[Abstract/Free Full Text].


Am J Physiol Endocrinol Metab 280(1):E179-E186
0193-1849/01 $5.00 Copyright © 2001 the American Physiological Society