## Abstract

To correctly evaluate the glucose control system, it is crucial to account for both insulin sensitivity and secretion. The disposition index (DI) is the most widely accepted method to do so. The original paradigm (hyperbolic law) consists of the multiplicative product of indices related to insulin sensitivity and secretion, but more recently, an alternative formula has been proposed with the exponent α (power function law). Traditionally, curve-fitting approaches have been used to evaluate the DI in a population: the algorithmic implementations often introduce some critical issues, such as the assumption that one of the two indices is error free or the effects of the log transformation on the measurement errors. In this work, we review the commonly used approaches and show that they provide biased estimates. Then we propose a novel nonlinear total least square (NLTLS) approach, which does not need to use the approximations built in the previously proposed alternatives, and show its superiority. All of the traditional fit procedures, including NLTLS, account only for uncertainty affecting insulin sensitivity and secretion indices when they are estimated from noisy data. Thus, they fail when part of the observed variability is due to inherent differences in DI values between individuals. To handle this inevitable source of variability, we propose a nonlinear mixed-effects approach that describes the DI using population hyperparameters such as the population typical values and covariance matrix. On simulated data, this novel technique is much more reliable than the curve-fitting approaches, and it proves robust even when no or small population variability is present in the DI values. Applying this new approach to the analysis of real IVGTT data suggests a value of α significantly smaller than 1, supporting the importance of testing the power function law as an alternative to the simpler hyperbolic law.

it is widely recognized that an accurate and exhaustive evaluation of glucose regulation in a subject must take into account both the ability of insulin to lower glucose level (insulin sensitivity) and the effect of glucose in promoting insulin secretion (β-cell responsivity). These two processes must be interpreted jointly, because it is only the balance of the two that can give an informative picture of the actual status of the system.

The disposition index (DI) is, to this day, one of the most widely used approaches to do so. Although other authors (25) had previously detected a similar relationship, the paternity of the original formulation is attributed to Bergman et al. (9). This paradigm assumes that subjects with similar efficiency of the glucose disposal metabolic system may have different values of insulin sensitivity [SI; e.g., provided by the intravenous glucose tolerance test (IVGTT) glucose minimal model (8)] and β-cell responsivity Φ, [e.g., provided by the C-peptide minimal model (24)] but share a similar value of their product,

DI is thus a mathematical tool to evaluate the efficiency of glucose homeostasis in a given individual (β-cell compensation), which is the ability of the β-cells to upregulate insulin secretion in response to a decrease in insulin sensitivity. Due to the implied inverse proportionality between SI and Φ, this is widely known as the hyperbolic law. More recently, a more comprehensive paradigm has been tested by Kahn et al. (15) with an exponential parameter α:

This parameter introduces an additional degree of freedom that is able to accommodate for a shape of the relationship between secretion and sensitivity indices (in the range of observed values) that is different from the hyperbolic law. An estimated value of α statistically different from one would allow rejecting the null hypothesis that the relationship is hyperbolic. Hereinafter, we will refer to the alternative DI paradigm of *Eq. 2* as power function law.

In the original study by Bergman et al. (9), the relationship between minimal model-derived insulin sensitivity and second-phase posthepatic insulin secretion was considered, but ever since then, many different indices (especially for insulin secretion) have been proposed for the calculation of the DI. A large number of studies in the literature (7, 13, 15, 16c, 17, 19, 26–28) employed the DI and addressed the problem of its determination. They also inspected, for different combinations of sensitivity and secretion indices, whether the parameter α is significantly different from 1, suggesting that the relationship may be more complicated than a rectangular hyperbola. Traditionally, the most widely used approach consists in considering insulin sensitivity and β-cell responsivity as coordinates of a Cartesian plane and then employing a traditional fit approach to find the curve (hyperbola or power function) best explaining the data. However, to avoid the difficulties of a nonlinear two-variable fit, simplifications such as the log transformation of the data or assumptions on the precision of the individual indices have been used in the literature. The first contribution of this work is a novel fit approach based on a nonlinear total least squares (NLTLS) technique that does not employ similar approximations. With the use of simulated data, this new algorithm is tested and compared with the other ones used previously in the literature, and, as shown in results, it provides much more accurate results.

Similarly to the all the other curve-fitting approaches, NLTLS implicitly relies on the hypothesis that all the subjects in the population have the same value of DI. Thus, the only variability that is assumed to be present in the data is the estimation uncertainty for the indices of insulin sensitivity and β-cell responsivity in each subject. This is arguably an oversimplification; even if a group of subjects is classified as healthy, it is very unlikely that they will all share exactly the same value of DI. As shown in this work, this assumption is extremely important, and its consequences should not be overlooked.

Therefore, the second contribution is a population approach based on a nonlinear mixed-effects model (NLMEM) that accounts separately for between-subject variability (BSV) in the DI values and the residual unexplained variability (RUV) in the data. The BSV quantifies the physiologically driven variability between the subjects in the study, whereas RUV (a term commonly used in the NLMEM jargon) refers to other sources of variability that cannot be ascribed to BSV, i.e., model error and measurement error, the latter being, in this case, the uncertainty affecting the estimates of individual indices of sensitivity and secretion. In the new framework we propose, the population-typical values of DI and α are considered as features that characterize the population distribution of the indices of insulin sensitivity and β-cell responsivity. The NLMEM analysis can be carried out with modeling tools such as NONMEM (see Ref. 4 for the *NONMEM Users Guide*) or SPK (23), which we employed here, but implementations in R (18), Monolix (16b), or WinBUGS (16) are also available.

All of the methods explained in the present work, including the two novel approaches (NLTLS and NLMEM), are compared on simulated data sets, which were generated according to two DI situations, the former with all subjects sharing the same DI value and the latter with variability between subjects in the DI values, but reproducing in both cases the variability observed in real data (2, 3). Finally, results on real data are also presented.

## METHODS

### Real and Simulated Data Sets

#### Real data.

The real data set derives from a previous study (2, 3) and consists of 204 healthy subjects who underwent an insulin-modified IVGTT with a full sampling schedule (240 min, 21 samples). The individuals belong to two large groups characterized by different ages: 59 young (age 23 ± 3) and 145 elderly subjects (age 69 ± 6). For each of them, identification by the modeling software SAAM II (1) of the IVGTT glucose (8) and C-peptide (24) minimal model provided individual estimates of both the value and precision of insulin sensitivity index (SI) and first- and second-phase secretion indices Φ_{1} and Φ_{2}, which were subsequently combined into the total index Φ_{tot}.

Both indices were approximately log-normally distributed across the patients in the population. Their precision was evaluated by the parameter estimation algorithm implemented in SAAM II by propagating the measurement errors of glucose and C-peptide concentrations. Results on our 204 subjects indicated that precision of SI and Φ estimates are well described by a constant SD and constant coefficient of variation (CV) model, respectively, and the variability of their precision across the population, in both cases, followed a log-normal distribution. The SDs for SI had an average value equivalent to 5% of its population typical value, whereas the CVs for Φ were on average ∼8%. The structure and size of the variability detected in the real data were the basis for the generation of the simulated data explained below.

#### Simulated data.

Since the true value of DI in our population is unknown, real data allow only a comparison between the results of the different methods. To assess their real performance, simulated SI-Φ data were generated in keeping with the statistical properties of the real data and satisfying predefined DI laws, namely the hyperbolic (*Eq. 1*) and power function laws (*Eq. 2*).

The simulation comprised two steps; first, the “true” values of SI and Φ were generated for each subject, and second, these “true” values were perturbed to mimic noisy experimental data. As regards the first step, data sets were generated either with no variability in DI values, i.e., true values of the indices satisfy either *Eq. 1* or *2*, or with variability between subjects in the DI values. The general model for the true values of SI_{i} and Φ_{i} is a joint population distribution (see appendix for the rationale):
^{2} a scaling parameter for the variance. The parameter α changes the shape of the DI curve, whereas ω^{2} regulates the size of SI and Φ between-subject variability.

As regards the second step, measurement errors to be added to the true values were generated by first simulating the precisions (constant SD for SI and constant CV for Φ) from a log-normal distribution and then using the selected precision as variability of a Gaussian distribution to generate the actual error (truncated to prevent negative values).

#### No population variability in the DI values.

A first data set was generated with no variability in the DI values. This was done by fixing the parameter ρ to −1, since, as shown in the appendix, a perfect anticorrelation between the two indices (ρ = −1) forces each individual deviation of SI from the mean value to be compensated entirely by an opposite deviation in Φ so that all subjects share the same DI value. The value of ω was fixed to ∼0.5 to reproduce the variability of ∼50% for both indices observed in the real data. The average values θ_{SI} and θ_{Φ} were tuned to obtain 100 as a symbolic expected value for DI. Three different values of the parameter α were used (0.5, 1, and 2), and for each of them, 100 sets consisting of 1,000 subjects were generated. For each set, the true DI values all superimpose exactly over the same DI curve. When perturbed according to the precision estimated in the real data, their level of scattering was visibly much smaller than in the real data set. A possible explanation could be that the estimated precision of SI and Φ reflects only measurement errors and not other sources of variability, e.g., model error and interday variability. Thus, even if the measurement noise was generated according to the above-mentioned constant CV and SD models, the level of uncertainty used in the simulation was increased by about five times (40% CV for Φ and SD equivalent to 25% of the typical value for SI) to produce simulated data in good agreement with the real ones, which were indeed extremely scattered.

#### Population variability in the DI values.

In a second data set, the assumption that all simulated subjects share the DI value was removed, and some BSV for DI_{i} was introduced. SI_{i} and Φ_{i} were still sampled from the joint distribution in *Eq. 3* but with progressively increasing levels of variability in the DI values, which were obtained by changing the value of parameter ρ from its initial value, −1 (−0.8, −0.5, and −0.2 were used). At the same time, the RUV was decreased, maintaining the overall variability at a similar level (with respect to the precision values estimated in the original data set for the secretion and sensitivity indices, the levels of uncertainty were multiplied by 4, 1, and 1/5, respectively). In this way, the focus of the analysis was on the consequences of the different hierarchy underlying the variability rather than its size. The BSV for SI and Φ (dependent only on the values of parameters ω and α) remained the same in these simulations, but due to the alteration of the value of the correlation ρ, the BSV in the DI changed; e.g., for the case α = 1 it was equal to 32 (ρ = −0.8), 50 (ρ = −0.5), and 63% (ρ = −0.2). Since BSV in the DI depends both on ρ and α, as shown in the appendix, the values are lower for larger α-values and vice versa.

### Parameter Estimation Methods

#### Linear fit on log-transformed data.

Curve-fitting approaches are based on the interpretation of the individual values of SI_{i} and Φ_{i} as coordinates *x* and *y* of a Cartesian system. They assume that the true (error-free) data lie on a curve, on a hyperbola, as in *Eq. 1*, or on a power function, as in *Eq. 2*, depending on the value of α. Finding the values of DI and α characterizing the population consists in identifying the curve that better fits the data, i.e., the curve lying at the minimum distance from the data points. The various methods differ on how they define this distance; many alternatives have been proposed in the literature, mostly employing approximations. Some approaches operate a log transformation to turn the problem into the fit of a straight line. Here, we will assess the most popular approach, which fits a straight line on log-transformed data with a one-variable approach; i.e., it assumes that only the *y* variable is affected by noise but uses a coefficient to adjust for the error in the second variable. The standard reference for this method is from Riggs et al. (20), but the first mention we could find is in a publication by Gini (14); thus it will be denoted as Log-Gini in this work. We will also test a linear total least squares approach applied after log transformation of the data (Log-TLS). Details on both methods are included in the appendix.

#### NLTLS.

The first novel method proposed here, named NLTLS, extends to the nonlinear case the minimization of the sum of squares of the weighted residuals along both directions, thus considering the presence of errors in both variables. To this purpose, for each data point (SI_{i}, Φ_{i}), its projection on the curve, denoted as (SI_{i}*, Φ_{i}*), is defined as the point lying on the curve and having the shortest distance to the data point, where the distance is calculated by weighting the components along *x* and *y* with the respective uncertainties (i.e., the estimated precisions σ_{SIi}^{2} and σ_{Φi}^{2}). Residuals in the TLS fit are equal to these distances; thus the objective function, which is minimized with respect to DI and α, can be written as follows:

As can be seen in *Eq. 4*, uncertainties along the *x*- and *y*-direction are properly taken into account by weighting the residuals along *x* and *y*, res_{x} and res_{y}, with the respective variances derived from the estimation uncertainties of the sensitivity and secretion indices. Efficient algorithms are available for the TLS fit of a straight line, but in the case of a power function we were not able to identify an analytical solution to the calculation of the projection (SI_{i}*, Φ_{i}*). Thus a numerical minimization procedure was applied, which yields the result at the price of a longer computation time. Within the minimization procedure for the objective function in *Eq. 4*, a further level of minimization was required to compute the projections (SI_{i}*, Φ_{i}*) of each data point (SI_{i}, Φ_{i}) on the curve. Matlab (16a) was used for this implementation, and code is available upon request.

#### NLMEM approach.

All the methods proposed previously, including NLTLS, are based on a curve-fitting approach, and therefore, they rely implicitly on the existence of a single curve (and thus DI value) that characterizes all subjects in the population. The variability in the data is assumed to be due only to the uncertainty with which the sensitivity and secretion indices are measured. This assumption makes these approaches inconsistent in the case of BSV in the DI values. Our newly proposed approach instead accounts for the presence of BSV in the observed DI values by postulating a joint probabilistic description of SI_{i} and Φ_{i} similar to *Eq. 3* and estimates the unknown parameters θ and Ω by using a population approach, such as NLMEM (5, 10). A model was implemented in NONMEM VI 2.0 and used to estimate the population parameters θ_{SI}, θ_{Φ}, ω, α, and ρ of *Eq. 3* and consequently combine them as detailed in the appendix to obtain the typical value of the DI. In the estimation process, the estimates of individual indices SI_{i} and Φ_{i} were used as data together with their precision. The NONMEM code for the implementation of this model is available upon request.

#### Comparison.

All methods were assessed on simulated data sets; for each setting of α and ρ used in the simulation, the estimates of the 100 repetitions were collected and displayed in boxplots, which show the median value and the size of the quartiles. The median can be used to assess the accuracy of the estimation algorithm, whereas the width of the boxes and whiskers is an indicator of its precision. For the NLTLS approach, we also collected the value of σ̂^{2} (the variance a posteriori of the weighted residuals), which was calculated with the following formula:
*n* is the number of data points (and thus *2n* is the number of residuals), and p is the number of parameters in the model, in this case 2, DI, and α.

## RESULTS

### No Population Variability in the DI

Results obtained on the data sets with no BSV in the DI values are summarized in Fig. 1. The estimates provided by the methods using a log transformation are affected by some bias, in particular the DI estimate. Log-Gini tends to underestimate DI (about −40% or worse), whereas Log-TLS overestimates it (about +40% or more). The precision and accuracy with which α is estimated are more satisfactory, and Log-TLS performs slightly better (−15% for Log-Gini, about +10% for Log-TLS). The newly proposed NLTLS and NLMEM methods yield the best performance in terms of both precision and accuracy and for both parameters. The size of the bias for DI is significantly lower compared with the other approaches, at most ∼20% as opposed to at least ∼40%. A similar trend is also evident for the precision of the estimates. In fact, NLTLS and NLMEM results are generally less scattered than the other algorithms. The improvement in the estimate of α is also significant; the bias is at most ∼5% instead of ∼10% for Log-TLS and ∼20% for Log-Gini. The parameter α used in the simulation (0.5, 1, and 2) does not seem to affect the outcome of the comparison.

### Population Variability in the DI

Results on data characterized by BSV in the DI values are shown in Figs. 2⇓–4. Figure 2 reports the results for a situation (*ρ* = −0.8) in which the variability in the data is still due mainly to the uncertainty in the estimation of the individual indices rather than the effect of BSV, which is limited. NLMEM provides slightly more precise estimates, but NLTLS and Log-TLS also yield satisfactory results. Interestingly, the performance of the latter two methods seems affected by the value of α, whereas this phenomenon is not observed for NLMEM. As the level of BSV in the DI increases, as depicted in Figs. 3 and 4 (ρ = −0.5 and ρ = −0.2, respectively), the estimates of α and DI provided by the curve-fitting approaches become less and less reliable. This is expected, since when a large fraction of the variability is caused by differences between subjects (BSV) rather than estimation uncertainty of the individual parameters (RUV), the hypothesis underlying the curve-fitting approaches becomes a gross simplification, resulting in a large bias in the estimates of both α and DI. It is very interesting to notice that the increasing contribution of BSV seems to affect mostly the accuracy of the estimates rather than their precision. The method that seems to be more affected by the introduction of BSV is Log-Gini, whose bias in the estimate of both DI and α is significantly larger than the other methods.

It is also interesting to consider the estimates of σ̂^{2} (the variance a posteriori of the weighted residuals) obtained with the different data sets and reported in Table 1. The value of σ̂^{2} can be used as a marker for goodness of fit. If the a priori knowledge of the level of error in the data is reasonable and the model is correct, a value of ∼1 is expected for σ̂^{2}. As can be seen in Table 1, the values are below or ∼1 for the data sets with no or relatively small population variability in the DI (ρ = −1 or −0.8). Instead, much larger values are obtained when the BSV of the DI is larger and the RUV is smaller (ρ = −0.5 or −0.2).

### The Real Data Set

Finally, both NLTLS and NLMEM were applied to the real data sets. Since NLTLS is not designed to deal with large variability in the population and the fit can become very sensitive to outlier subjects, a jackknife technique (12) was employed to detect and remove such values. On average, the data of one or two subjects were removed from each data set. The analysis was first executed on the young and elderly populations separately and finally on the whole group. In this last case, BSV in the DI is present by hypothesis (since young and elderly patients are expected to have different DI values), so the NLTLS approach was used only as a mean of comparison, whereas the results are somehow expected to be unreliable. Table 2 reports the parameter estimates and their 90% confidence interval provided by a 3,000-sample bootstrap employing the NLTLS and NLMEM algorithms for DI and α and for SI coupled with Φ_{tot} as well as with first- and second-phase responsivity Φ_{1} and Φ_{2}. Not surprisingly, some of the results differ significantly between the two methods for both DI and α. However, the validity of the NLTLS estimates is questionable due to the values of σ̂^{2} detected by NLTLS, which are always much larger than 1, suggesting that the amount of unexplained variability in the data is much greater than indicated by the precision of the individual estimates. This might suggest that either the estimation confidence intervals for the individual indices SI and Φ are heavily underestimated or the model is inadequate to describe the data, or there is another source of variability that is not properly accounted for (e.g., BSV). This is also evident in Figs. 5, 6, and 7, which provide a graphical representation of the data and the NLMEM fits for the separate groups of young and elderly subjects. The size of the error bars (representing the RUV) does not seem large enough to explain the distance of the individual values from the DI curve. A further confirmation of this comes from the results provided by NLMEM, which detects values of correlation *ρ* that are consistently different from −1 (about −0.5), indicating a remarkable amount of scattering of the DI values, particularly when the group of elderly subjects or the whole population with the first-phase index Φ_{1} as β-cell responsivity index is analyzed. Estimates of power function law parameter α are <1 for all the data sets, but the difference is statistically significant only for the elderly cluster or when the entire population is employed. However, it should be pointed out that the cluster of young subjects is smaller (59 vs. 145), and therefore, it is harder to achieve statistical significance.

In a NLMEM framework, another statistic that is used to assess the significance of an additional parameter in a model is the drop in objective function value, which is approximately *χ*-square distributed and can therefore be used to obtain a *P* value. The *P* values for the additional parameter α contained in Table 2 confirm the same trend outlined by the confidence interval; the exponent is significantly smaller than 1 for the elderly cluster or when young and elderly subjects are analyzed concomitantly. Whichever method is used to assess the statistical significance, a clearly different pattern emerges in elderly vs. young subjects; DI decreases, suggesting an impairment in β-cell ability to compensate a reduction in insulin sensitivity, but also α decreases, suggesting that the β-cell impairment is more evident in individuals with lower SI.

## DISCUSSION

In this study, the performance of conventional and novel approaches to evaluate the DI law in a population has been assessed on both real and simulated data sets. When no population variability in the DI values is used for the simulations, the novel approaches NLTLS and NLMEM outperform the more conventional log transformation methods.

The likely reason is that although the use of a log transformation radically simplifies the problem, it introduces a nonlinear rescaling of the axes and thus of the residuals, i.e., the distances between each data point and their model predictions, and this can substantially change the leverage that each point has in the fit procedure. Moreover, as already mentioned, the log transformation also affects the error structure; e.g., a constant CV structure in the original data results in a constant SD of the log-transformed data. It is hard to tell which of the two phenomena affects the results more severely; the nonlinear rescaling of the data is expected to underestimate DI, since the subjects with high SI or Φ (and DI) values tend to cluster together, whereas the distance between the points close to the axes (low SI or Φ and low DI) increases and drives them farther from the average curve.

An example to clarify this phenomenon can be found in the difference between the geometric and arithmetic mean of a set of positive numbers. The calculation of the geometric mean is equivalent to calculating the arithmetic mean after log-transforming the numbers. And it is known that, given any set of positive numbers (denoted as z_{i}) that are different from each other, the geometric mean is lower than the arithmetic mean:

On the other hand, the first-order approximation used for the error structure is expected to have the opposite effect, partially compensating for the uneven stretching. Using the CVs of the original values as SDs in the log-transformed space (where the fit is performed) gives a relatively stronger leverage to the high values of SI and Φ. Possibly, this is why Log-Gini, which assumes the same degree of precision for all of the data points in the log-transformed plane and thus is not prone to this latter effect, tends to underestimate DI. Log-TLS, on the other hand, tends to overestimate DI, whereas it provides a rather accurate estimate of α.

NLTLS is computationally more cumbersome but provides reliable results because it does not use approximations to calculate the distance between the data points and the curve. NLMEM also performs well, even if it represents an overparameterization of the problem, in the no-population variability case.

The situation is different if population variability is present in the DI values; as already mentioned, the presence of other sources of variability in the data besides the estimation uncertainty of the secretion and sensitivity indices undermines the validity of all of the methods considered here, except NLMEM. And indeed, the performance of all of the curve-fitting approaches degrades significantly when population variability in the DI values comes into play, even if the overall amount of variability in the data is not altered. The underlying reason is that, as the magnitude of population differences (BSV) surpasses the effect of individual uncertainty (RUV), the estimates of precision for SI_{i} and Φ_{i}, even if correct, become less and less informative about the actual distance of the data points from the curve characterizing the typical DI of the population. In fact, the precision on the individual parameter estimates (σ_{SIi}^{2} and σ_{Φi}^{2}) does not carry any information about how scattered the DI values are in the population. A subject may be characterized by very precise SI and Φ estimates, but the resulting DI may be extremely different from the population DI typical value. The curve-fitting approaches are not designed to account for this kind of hierarchical variability, so they fail when different sources of variability come into play, even if no approximation is used, as in the case of NLTLS. A symptom of this is the value of σ̂^{2}; values much larger than 1 are detected when BSV is used to simulate the data. This shows lack of fit and indicates possibly that the a priori level of uncertainty is not large enough to account for all the variability detected a posteriori in the data. As a consequence of poor fit, an extremely large bias is introduced in the parameter estimates, whereas the level of precision is similar. Moreover, the size and direction of the bias becomes strongly dependent on the value of α used in the simulation of the data set. The different settings of α introduce asymmetries in the data that lead to systematic over- or underestimation. On the contrary, NLMEM is designed to capture the hierarchical structure of the variability and provides reliable results regardless of the value of α.

Finally, we have some comments on the results obtained using NLTLS and NLMEM on the real data (Table 2). The NLTLS fit approach detected a large value of σ̂^{2}, and NLMEM estimates levels of correlation ρ that are different from −1 among the individual indices. Both of these markers are indicators of the presence of considerable variability not explained in terms of precision of the individual estimates, suggesting that NLTLS results should be regarded as unreliable. Therefore, with the purpose of investigating the statistical significance of parameter α and thus the validity of the power function vs. the pure hyperbolic law in our data set, the NLMEM results are used.

The NLMEM results suggest that the typical DI values for the elderly population are lower than for the young one, and the graphical representations in Figs. 5⇑–7 show that the curve for the elderly subjects is closer to the origin than the one for the young ones.

Moreover, although the analysis of other data sets would be needed to confirm the current findings, the estimates of α provided by NLMEM are often lower than 1, supporting the importance of testing the power function law as an alternative to the simple hyperbolic law. A value of α less than 1 produces a curve rotated anticlockwise, with values lower than a rectangular hyperbola for low values of SI, thus suggesting that compensatory mechanisms are less efficient in subjects characterized by low insulin sensitivity.

### Conclusions

In this work, we showed the importance of the regression methodology in the analysis of the DI. We pointed out the shortcomings of the currently used methods and proposed the novel NLTLS fit. Although this new approach proved more reliable than its competitors, our results highlighted the inadequacy of all curve-fitting approaches in the case of variability in the DI values across the population. Analysis of real data confirmed that a considerable amount of BSV in the DI values is probably present, even in a population of healthy individuals. Thus we designed and tested a novel population model that was able to efficiently account for BSV. We recommend the use of a population approach in the study of the DI, since neglecting to account for the hierarchical structure of the variability can yield very unreliable results.

In our analysis, a two-step approach was used. First, the individual values of insulin sensitivity and secretion were calculated from the raw blood sampling data, and then the DI analysis on the individual results was performed. A more complete approach would consist in integrating the two steps into a joint population model. Nonlinear mixed-effects modeling has already been successfully applied with the glucose minimal model to estimate SI (11), and similarly it could be used for the estimation for insulin secretion. Thus, one can imagine applying such an approach to jointly estimate the secretion and sensitivity indices directly from glucose and insulin profiles and, at the same time, extracting the information on DI and α from the population covariance matrix.

## GRANTS

P. Denti was partially funded by the scholarship “Borsa Gini”, sponsored by the Gini Foundation, Padua, Italy.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the authors.

## AUTHOR CONTRIBUTIONS

P.D. analyzed the data; P.D., G.M.T., and C.C. interpreted the results of the experiments; P.D. prepared the figures; P.D. drafted the manuscript; P.D., G.M.T., and C.C. edited and revised the manuscript; P.D., G.M.T., and C.C. approved the final version of the manuscript; G.M.T. and C.C. did the conception and design of the research.

## ACKNOWLEDGMENTS

These results were presented partially at the EASD 2007 Annual Conference in Amsterdam, The Netherlands, and included in the conference proceedings (Denti P, Campioni M, Toffolo GM, Rizza R, Cobelli C. “Is the disposition index law hyperbolic? Importance of the regression methodology.” Abstract 602) and at PAGE 2009 in St. Petersburg, Russia (Denti P, Salinger D, Vicini P, Toffolo GM, Cobelli C. “A NonLinear Mixed-Effects Approach to the Estimation of the Glucose Disposition Index.” PAGE Abstract 1648).

P. Denti thanks Eng. Domenico Dei Tos for performing the analysis on the real data set as part of his Laurea Thesis project, Dr. David Salinger for advice on the NONMEM implementation of the population model, and Dr. Bradley Bell and Paolo Vicini for fruitful discussion.

Present address of P. Denti: Division of Clinical Pharmacology, University of Cape Town, Cape Town, South Africa.

## APPENDIX

#### Formulation of the DI as a Joint Population Distribution

If SI and Φ are assumed to be log-normally distributed, their joint population distribution can be written in the following general form:
_{i} and Φ_{i} and Ω is their covariance matrix expressed, without loss of generality, in terms of variances (ω^{2}) and coefficient of correlation (ρ). Under the assumptions that both SI and Φ are log-normally distributed, the distribution of the DI is also log-normal, with mean and variance related to mean values and variances of the sensitivity and secretion indices, as specified in the following equation:

From *Eq. 8*, the conditions under which all subjects share the same DI value, i.e., the variance of DI is equal to zero, are as follows:

Apart from the trivial solution ω_{SI} = ω_{Φ} = 0 (i.e., sensitivity and secretion indices assume the same value in all subjects), *Eq. 8* is satisfied if
_{i}) and log(Φ_{i}) are perfectly negatively correlated. Moreover, the ratio between their population variances determines the value of the parameter α.

Therefore, to interpret the SI-Φ data with the DI paradigm, we can use the following easy reparameterization

The value of the correlation ρ can then be used to determine the size of BSV in the DI values.

#### Log Transformation

The most widely used curve-fitting approach, introduced by Kahn et al. (15), consists in applying a log transformation to the original data, turning the problem into a linear fit. In fact, the log transformation of both variables, SI and Φ, reshapes the hyperbola into a straight line, as shown in the following formulae:

However, this transformation is nonlinear and does not preserve the Gaussian error structure of the model; in this way, the hypotheses supporting the consistency of the classical fit approach are not satisfied. To maintain at least part of the information on the precision of the individual estimates, a first-order approximation can be used. Loosely speaking, this is equivalent to assuming a Gaussian distribution of the error for the log-transformed estimates, using the values of the CV of the original estimates. However, this simplification is quite poor, especially if the uncertainty is large, since the nonapproximated probability distribution of the log-transformed data is nonsymmetric and skewed to the right. On top of this, one should bear in mind that there is no guarantee that the optimal solution in the log-transformed space preserves its optimality under nonlinear transformations. Therefore, especially if a large variability is present, the two solutions might be very different.

In addition, even if the log transformation simplifies the problem to the fit of a straight line, the error is still present in both variables, and a proper technique must be employed to take this into account. Many approaches are used in literature, but the most widespread is the perpendicular weighted regression presented by Riggs et al. (20), which is not, however, the original source. The oldest reference to the same formula we found is in Gini (14), so we denoted it Log-Gini. This method is essentially an orthogonal fit that uses a correction factor λ (ratio of the variances) to account for the difference of precision along the two directions. Therefore, technically, a fit along just one direction is performed, but the correction factor λ allows one to also consider the second direction. Although a fast linear method is available, and thus no minimization of the objective function is explicitly required, we propose here the formula for the sake of comparison with the other methods:

The symbol ⊥ denotes the orthogonal projection on the straight line once the axes have been rescaled to accommodate for the heteroscedasticity along the different directions. However, a limitation of Log-Gini is represented by the fact that it assumes the same degree of precision for all the subjects. In this way, all the data points will be considered equally reliable in the fit. Since in our simulations we included individual precisions, a median value was employed to obtain an estimate of λ, as shown in the following formula:

An alternative approach to execute the fit for the straight line after the log transformation has been performed is TLS:
_{i}*), log(Φ_{i}*)) is the prediction in the TLS sense, i.e., the point on the straight line lying at the shortest weighted distance from (log(SI_{i}), log(Φ_{i})). As mentioned already, in this linear case, the TLS projection can be obtained with a closed-form formula, and therefore, it is computationally very fast. This approach is referred to here as Log-TLS.

- Copyright © 2012 the American Physiological Society