## Abstract

The minimal model method is widely used to estimate glucose effectiveness (S_{G}) and insulin sensitivity (S_{I}) 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 β-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 S_{G} and S_{I} 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 S_{G} = 2.66 ± 1.09 × 10^{−2} · min^{−1} and S_{I} = 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 S_{G} and S_{I}. However, RSS and ITS gave S_{G} = 2.59 ± 0.73 and S_{I} = 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

the minimal model of glucose disappearance (5) is widely used to estimate metabolic indexes of glucose effectiveness (S_{G}) and insulin sensitivity (S_{I}) 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 S_{I} 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 (S_{I}) 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 S_{G} and S_{I} 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 S_{G} and S_{I}. 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 β-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

#### 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
Equation 1
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, G_{b} (= Q_{b}/V) and I_{b} 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: S_{G} (min^{−1}), S_{I}(min^{−1} · μU^{−1} · ml), p_{2} (min^{−1}), the insulin action parameter, and V (dl kg^{−1}), the glucose distribution volume. S_{G} and S_{I} 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, S_{G} measures the ability of glucose per se, at basal insulin, to stimulate glucose disposal and to inhibit glucose production. Similarly, S_{I} 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** = [S_{G}, S_{I}, p_{2}, 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 S_{G} and S_{I} consists of estimating S_{G}and S_{I} separately in each subject via nonlinear regression and then calculating the sample mean (μ) and covariance (Σ^{2}) of all the S_{G} and S_{I}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 Σ^{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, **p**
_{j}
Equation 2where N_{j} is the number of data points available for the j^{th} subject, t_{i,j} and G_{i,j}
^{OBS} are the i^{th} time point and data point, respectively, of the j^{th} subject, ς^{2}
_{i,j} is the variance of the measurement error of the i^{th} data point, and G(**p**
_{j},t_{i,j}) is the minimal model prediction of glucose concentration for a given**p**
_{j}.

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), V_{j}, of the resulting estimate **p̂**
_{j} of **p**
_{j}from the inverse of the appropriate Fisher Information Matrix
Equation 3where the superscript T indicates vector or matrix transpose, and [∂G/∂**p**
_{j}] 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 (G_{b} and I_{b}) as the end-test basal concentrations, i.e., the last available measurement (usually at 240 min). Figure1 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).

We calculated the population mean for each parameter as the sample mean of all the individual parameter estimates (**p̂**
_{j})
Equation 4where N is the number of subjects, and we calculated the population variance as the corresponding sample variance
Equation 5From a methodological standpoint, this approach is undesirable for a number of reasons. It neglects the fact that the precision of the estimates of S_{G} and S_{I} 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 V_{j}, 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}, Σ^{2}(0) = Σ_{STS}
^{2}.

*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 **p**
_{j}
Equation 6where the distance of the current parameter estimate from the population mean (the prior) is also penalized; we denote with**p**
_{j,i} the i^{th} element of the parameter vector for subject j and with N_{p} the number of elements in the vector **p**
_{j}. The estimate**p̂**
_{j} obtained by minimizing this objective function is often called post hoc, or empirical Bayes, estimate. We can again calculate V_{j}, the precision of the parameter estimate**p̂**
_{j}, from the Fisher Information Matrix (see below); μ_{i}(k) is the value of the population mean at the k^{th} iteration of the method, and Σ_{i,i}(k) is the i^{th} diagonal element of the population covariance matrix at the k^{th} iteration. Calculate, then, the updated population mean of the parameter vector
Equation 7and the covariance
Equation 8
*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) Equation 9

where Σ^{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 Σ^{2} is a positive definite matrix by construction), the question arises of how to choose a proper value for Σ^{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 S_{I} and S_{G}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

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.

We first applied the STS method to both FSS and RSS data. We detected a significant difference (*P* < 0.05) between S_{I} 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 S_{I} is fairly stable (except when its value is very low, as in *subject 11*), whereas S_{G} 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.

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 S_{G} and S_{I} was always between 0.97 and 1.00.

## DISCUSSION

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 S_{G} and S_{I}(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 S_{I} 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 Σ for p_{2} 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 ς^{2}
_{i,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)
Equation 10
*Equation 10
* gives twice the negative logarithm of the likelihood function of the data. If ς^{2}(p_{j},t_{ij}) 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 S_{G} and S_{I}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.

## Acknowledgments

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.

- Copyright © 2001 the American Physiological Society