INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION REFERENCES

THE MINIMAL MODEL OF GLUCOSE KINETICS, in conjunction with the insulin-modified intravenous glucose tolerance test (IM-IVGTT), is widely used to measure insulin sensitivity in subjects with impaired glucose tolerance and type 2 diabetes in both clinical and epidemiological studies (1, 14, 17, 23). Although some of its assumptions, in particular the single compartment approximation to describe glucose kinetics, have been recently reexamined (7, 8), it is important to understand the performance of the classical minimal model when it is identified with the most advanced estimation techniques. Virtually all minimal model identification strategies employed in the literature resort to nonlinear least squares (NLS) estimation. NLS provides a point estimate of each model parameter and, by means of the Fisher information matrix, a measure of its uncertainty expressed as standard deviation (SD) can be obtained, from which a confidence interval can be determined. However, this approach has difficulties in handling possible asymmetries in the probability distribution of the estimates. In providing the estimates, NLS exploits only the experimental data (e.g., plasma glucose concentration time series in our case) and the knowledge of the measurement error statistics. NLS is thus a purely data-driven approach. One reported problem with this approach in minimal model studies is that, in a number of subjects, especially those with type 2 diabetes, insulin sensitivity (S_I) is calculated as S_I=0 with the minimal model program MINMOD (19). For instance, Saad et al. (23) found that S_I=0 occurred in 12 of 24 type 2 diabetic subjects, i.e., a prevalence of 50%. As a consequence, in population studies, histograms with a likely artificial bimodal pattern (i.e., a peak at S_I=0 and another peak at a positive S_I value) are obtained. At the same time, with more general purpose software packages, e.g., SAAM II or ADAPT (2, 10), S_I is estimated to be nonzero, and thus physically sound, but very small and with a very large uncertainty. As a result, negative values are included in the confidence interval. The above interpretative difficulties, often referred to as the "zero" S_I problem, call for the adoption of parameter estimation techniques more sophisticated than NLS.

The so-called Bayes approach is a methodology that can be employed to attack the parameter estimation problem as an alternative to NLS, but it is less adopted in practice. Bayesian estimation techniques exploit not only the experimental data but also the a priori information (also denoted in the following by the term "prior," as is usual in the statistics literature) on the unknown parameters of the model (27). In physiological modeling, this a priori information, e.g., nonnegativity and range of variability, can typically be made available from physical considerations and population studies (see Refs. 8 and 24 for two recent applications of Bayesian estimation in metabolism). The key advantage of Bayes estimation over NLS is that the former returns the entire a posteriori probability distribution of the model parameters, from which the marginal probability distribution of each parameter can be obtained. From this probabilistic quantity, a point estimate and its confidence interval can be computed. Notably, with use of Bayes estimation, confidence intervals are in general allowed to be asymmetric with respect to the point parameter estimates, in contrast to NLS, where the implicit assumption of parameter estimates to be Gaussian leads to confidence intervals symmetric by construction. The limited popularity of Bayes estimation in physiological modeling is likely due to its heavy computational burden, because often the a posteriori probability distribution of the model parameters must be numerically computed by the Markov chain Monte Carlo (MCMC) simulation approach (13).

In this article we develop a Bayesian parameter estimation strategy to identify the minimal model of glucose kinetics in diabetic subjects and implement it by MCMC. First, we only supply to the estimator a rough prior, i.e., the unknown model parameters are positive. Then, because more refined a priori information on S_I can be made available from the literature, we reidentify the model and show improved results. The database consists of IM-IVGTT performed in ten subjects. Whereas NLS provides S_I estimates and confidence intervals that are not acceptable in four cases, Bayes estimation always determines a nonzero estimate of S_I and, when resorting to the refined prior model, a credible confidence interval. The Bayesian approach we propose should prove useful to reanalyze large databases of epidemiological studies, such as GENNID (20), the Finland-US investigation of NIDDM genetics (FUSION) (26), and the Insulin Resistance Atherosclerosis Study (IRAS) (17), where the numbers of subjects are in the range of thousands.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION REFERENCES

Subjects, Protocol, and Data

Minimal Model of Glucose Kinetics

(1)

(2)

1

1

1

1

View larger version (16K): [in this window] [in a new window]   Fig. 1.   Cold minimal model of glucose disappearance. NHGB, net hepatic glucose balance; D(t), cold glucose dose; Ib, basal insulin dose; G(t), glucose concentration at time t; Q(t), glucose mass at time t; Up(t), peripheral glucose disposal at time t; ki, rate constants characterizing material fluxes (solid lines) or control actions (dashed lines).

NLS Identification

(3)

As in Ref. 1, we estimate, in each of the 10 subjects of our database, the unknown parameter vector  from the data vector y = [y₁, y₂,..., y_N]^T by weighted NLS, i.e.

(4)

where the weights w_i are optimally chosen, i.e., equal to the known measurement error SD (6, 9). The precision of each estimate, expressed by its SD, is obtained from the square root of the diagonal elements of the inverse of the Fisher information matrix.

Bayesian Identification

 y

The a posteriori probability function of  contains extremely insightful information. For instance, a point estimate, called the posterior mean of the unknown parameter vector, can be determined as the expected value of the vector |y

(5)

The a posteriori probability function of  also allows the derivation of the 95% confidence intervals [_L,_H], e.g., the interval between the 0.025 and the 0.975 quantile of the a posteriori probability distribution of .

The a posteriori probability density function of the parameter vector  can be obtained from the Bayes theorem as

(6)

In the right-hand side of Eq. 6, f_y(y) does not depend on  and, for us, is just a scale factor we are not interested in (13), f() is given (see the next paragraph), and f_y|(y|) (often referred to as the likelihood function of the data) can be obtained from Eq. 3 by exploiting the knowledge of measurement error statistics. For instance, when, as in our case study, v = [v₁, v₂, ..., v_N]^T is Gaussian with zero mean and covariance matrix _v, the random vector y| is Gaussian with covariance matrix _v and mean given by [h(t₁,), h(t₂,), ..., h(t_N,)]^T.

A nonnegativity prior on the four parameters. Prior information is a key additional ingredient of Bayes estimation, i.e., one needs to specify the a priori probability density function f() of Eq. 6, which summarizes our beliefs on model parameters "before having seen the data." From Eq. 6 one can note that the Bayesian approach coincides with maximum likelihood, if a uniform distribution from and + is assumed as a prior for all of the unknown parameters. In our minimal model studies, we know in advance that all the parameters (S_G, p₂, G_o, and S_I) are intrinsically nonnegative. Therefore, it is reasonable to choose a prior f() that excludes negative parameter values. For S_G, p₂, and G_o, nonnegativity is the only information embedded into the prior distribution, i.e., an expected range of variability is not declared. In formal terms, we state that each value of S_G, p₂, and G_o is a priori equally probable in [0, a] with a. The same poorly informative prior could obviously also be applied to S_I. One can note from Eq. 6 that, in this way, the posterior density function coincides, except for a constant factor, with the likelihood function in the domain of interest. As will be illustrated in RESULTS, this prior, albeit rough, allows a solution of the "zero" S_I problem, because a confidence interval, including only nonnegative values, is obtained. However, because one additionally experienced problem with NLS is that the confidence interval of S_I tends to be large, it is worthwhile trying to incorporate more knowledge into the S_I prior, as we will explain.

A refined S_I prior. Ideally, to define the a priori probability density function of S_I, one would need a large database, e.g., containing thousands of subjects, to extract it directly. Unfortunately, such a database is not available to us. However, several insulin sensitivity NLS studies in diabetic subjects are available in the literature that allow one to obtain an approximation of the actual distribution of the true S_I values. From these investigations, one knows that it is unlikely that S_I will exceed certain values. For instance, in Ref. 18, 30 type 2 diabetic subjects were studied, and the average S_I was 0.7 (10⁴ min¹/µU ml¹) with a standard error of 0.1. Of note is that, in this data set, no "zero" S_I problems were present. Thus one expects that S_I values greater than a certain threshold are less and less probable. At the same time, because S_I can be very small, although we do not know how small it can be in a certain subject before seeing his or her data, it is wise not to incorporate in the prior any beliefs about its variability at low levels. Thus we have assumed for S_I an a priori probability density function in which S_I values <2 × 10⁴ (min¹/µUml¹) are equally probable, while S_I values >2 × 10⁴ are less and less probable, according to a decreasing exponential law (with exponent equal to 10⁴ µU ml¹/min¹). The chosen threshold value of 2 × 10⁴ min¹/µU ml¹, even if not directly obtained from a distribution of S_I estimates in NIDDM subjects, is realistic and is supported by the literature, where it is unlikely to find insulin sensitivity values in NIDDM subjects >2 × 10⁴ min¹/µU ml¹.

The MCMC method. The goal is to obtain f_|y(|y), i.e., the joint probability density function of the unknown minimal model parameters, given the glucose data. This task is analytically intractable, both because of the complex relationships between parameters and data and because of the shape of the probability distributions involved. A solution suggested in the literature is to derive f_|y(|y) in sampled form by adopting a simulation strategy known as Markov chain Monte Carlo (13). MCMC methods are based on two steps. In the first step, a suitable Markov chain that converges (in distribution) to a target distribution, f_|y(|y) in our case, is generated. Then, a Monte Carlo integration step is performed to numerically compute the integrals of interest, e.g., Eq. 5. There are many MCMC methods in the literature; see Ref. 13 for a review. They differ from each other in the way the Markov chain is generated. However, all of the different strategies proposed in the literature can be traced back to the Metropolis-Hastings algorithm (15). In this work, we generated the Markov chain by a symmetric transition kernel that extracts a sample of the chain around the previous one (this scheme is usually referred to as "random-walk Metropolis") by use of independent uniform probability densities. The variance of such densities was chosen case by case. As far as implementation issues of the MCMC method are concerned, we note that convergence of the chain was assessed by the Raftery criterion (21), which is also known in the literature as binary-control. To compute the number of iterations necessary to estimate a posterior quantile from a single run of a Markov chain, a pilot analysis of output samples was used to fit a two-state Markov chain model, and from it one can calculate the length of the burn-in period and, then, the number of further iterations required to estimate quantiles of interest with the required accuracy. In the sequel, we have been required to estimate quantiles 0.025, 0.25, 0.5, 0.75, 0.975 with precision factors of 0.02, 0.05, 0.05, 0.05, and 0.02, respectively, with a probability of 0.95. From these quantiles, numerically robust confidence intervals can be calculated.

	RESULTS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION REFERENCES

Data

View larger version (37K): [in this window] [in a new window]   Fig. 2.   Glucose (units at left, interpolated with solid lines) and insulin (units at right, interpolated with dashed lines) plasma concentrations of the 10 subjects.

NLS Identification

View larger version (7K): [in this window] [in a new window]   Fig. 3.   Nonlinear least squares (NLS) and Bayes insulin sensitivity (SI) estimates in 10 type 2 diabetic subjects. Vertical line, point estimates; horizontal bars, 95% confidence intervals; NLS estimates were assumed Gaussian to determine their confidence interval as ±1.96 SD; for Bayes estimation, point estimates and confidence intervals were obtained from information depicted in Fig. 5 (see RESULTS, A refined SI prior).

Note that, because the covariance matrix of the measurement error is assumed to be known from the data, NLS coincides with maximum likelihood (ML). If it were assumed to be dependent on the model, then NLS and ML are no longer equivalent on a theoretical basis. However, in this case, the two results are virtually superimposable (unpublished data).

Bayesian Identification

A nonnegativity prior on the four parameters. The a posteriori probability density function of S_I estimated in each subject is reported in Fig. 4. As somewhat expected, in subjects 1, 2, 4, 5, 8, and 10, where NLS was successful, the posterior provided (in sampled form) by MCMC is concentrated around the NLS estimate. However, the new method reveals asymmetries in the distribution of the estimates, particularly in subjects 1 and 4. This is not surprising, because the symmetry of NLS confidence intervals reflects only the Gaussian assumption implicitly made on the estimates. Moreover, in subjects 3, 6, 7, and 9, MCMC analysis reveals that the marginal posterior densities have a high peak located at low and realistic S_I values but that they have a long tail. The confidence intervals, however, would be too pessimistic, because it is known from the literature that such high S_I values are unrealistic, e.g., it is impossible that S_I is, say, equal to 100 × 10⁴ min¹/µU ml¹ (see subject 3 in Fig. 4). So, to obtain a more credible confidence interval, it is necessary to adopt a more realistic S_I prior.

View larger version (20K): [in this window] [in a new window]   Fig. 4.   The a posteriori probability density function of SI for the 10 subjects (area under each curve is unitary) from Bayesian identification, with a nonnegativity prior on the 4 model parameters. On the x-axis, SI is measured in 104 min1/µU ml1. The y-axis quantity is unitless.

A refined S_I prior. Results for S_G, p₂, G₀, and S_I are reported in Table 2, together with the 95% confidence interval. Results for S_I are also displayed in Fig. 5, where the estimated a posteriori probability density function of S_I in each subject is reported (the axes are now different from those used in Fig. 4), and in Fig. 3, where a summary of Fig. 5 results is depicted to make the comparison with NLS easier. One can note, by comparing Figs. 5 and 4, how the threshold influences the tail of the posterior distribution of some subjects. One can also note that the posterior probability density function is zero in an interval at the right of S_I=0. This allows us to exclude, on a probabilistic basis, that S_I is greater than a chosen threshold. From the a posteriori density function of S_I, one can obtain the interval where the true value of S_I falls with 95% of confidence as the interval between the quantiles 0.025 and 0.975. This interval is displayed in Fig. 3, together with the posterior mean of S_I. The novel features of the approach we propose emerges in subjects 3, 6, 7, and 9, where NLS results were largely unsatisfactory with noncredible confidence intervals. One can see that the marginal a posteriori probability density function of S_I does not collapse to zero (this also applies to Fig. 4, but it cannot be visually appreciated due to the larger range of S_I values shown in the x-axis). In all of these subjects, the posterior mean of S_I is higher than the NLS estimate (markedly in subjects 3 and 7), because the minimum variance estimate takes into account the entire shape of the posterior and, so, also the presence in some subjects of the tails. More important, in each subject the confidence interval now lies entirely in the positive axis and allows the investigator to infer that the probability that the true S_I lies at the left of the displayed bar is 2.5%. In summary, a nonzero S_I value is detected with high probability by Bayesian identification in these very insulin-resistant individuals. This result is not always possible with NLS, where, given the uncertainty of the S_I estimate, one can only conclude that S_I is indistinguishable from zero (Fig. 3, e.g., subject 7). A credible confidence interval is obtained, which allows the minimal model user to determine with 97.5% of confidence a (nonzero) lower-bound number for S_I.

View larger version (20K): [in this window] [in a new window]   Fig. 5.   The a posteriori probability density function of SI for the 10 subjects (area under each curve is unitary) from Bayesian identification, with a refined SI prior. On the x-axis, SI is measured in 104 min1/µU ml1. The y-axis quantity is unitless.

	DISCUSSION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION REFERENCES

In this article, we have attacked the identification of the minimal model of glucose kinetics in type 2 diabetic subjects in a Bayesian context. The motivation for this stems from the interpretative difficulties generated from what is usually referred to as the "zero S_I" problem, i.e., in a significant number of cases NLS estimation returns an S_I=0 when the minimal model program MINMOD (19) is used, or very small values with an unrealistic confidence interval including negative values with software packages like SAAM II (2) or ADAPT (10).

In this case study, the Bayesian approach allowed the derivation (albeit only numerically through MCMC) of the joint a posteriori probability distribution of the model parameters, and thus the determination of a nonnegative point estimate of S_I, together with a plausible confidence interval. The differences of the point estimates between the Bayesian approach and NLS, at least in most subjects, are not large, but it is when the confidence intervals of the parameter estimates are considered that the Bayesian approach emerges as superior. In fact, because we think that the chosen prior as well as the confidence intervals employed are both reasonable, one can exclude on a probabilistic basis that S_I in the critical subjects is higher than a given threshold. Moreover, from the confidence interval obtained by Bayes estimation, one can exclude with 97.5% of probability that the true S_I lies at the left of the displayed variability range. In summary, the Bayesian approach is able to cope with those situations in which the confidence interval, obtained as two times the SD of the error given by the Fisher information matrix, is not realistic, because such a method is unable to cope with asymmetries in the distribution of the estimates, like those present in this study.

A critical point in Bayes estimation is the definition with a probability density function of the information available on the unknown model parameters before seeing the data. As far as S_G, p₂, and G₀ are concerned, we have limited ourselves to express the poorly informative knowledge that they are nonnegative, i.e., a uniform distribution along the positive axis was chosen. A more sophisticated prior was chosen for S_I, the very critical parameter of the model, based on the observation that, for values of S_I greater than a certain threshold, the higher the candidate S_I, the lower, in all likelihood, its probability. We have translated these beliefs about S_I in statistical terms by describing its a priori probability density function as a uniform distribution between zero and a certain threshold (i.e., 2 × 10⁴ min¹/µU ml¹, a value that appears as prudential on the basis of published studies on insulin resistance in type 2 diabetics) and an exponential decay for S_I greater than this threshold. It is important to note that, in the prior chosen for S_I, we do not include beliefs on how small it can be, i.e., where the whole "zero S_I" problem began. Choosing also for S_I a uniform distribution prior, as in Bayesian identification with only the nonnegativity prior, one would have obtained an S_I a posteriori probability density function lower than the one presented. This means that low S_I values are less probable and supports the robustness in each subject of the 0.025 quantile of Fig. 1, which allows us to exclude with 97.5% of probability that the true S_I lies at the left of the displayed range. Another consequence of using a uniform prior for S_I along the positive axis is represented by a gratuitous shift to the right of the 0.975 quantile. In fact, the new technique permits us to show that, in some cases, the model defines likelihood functions which assign a nonnegligible probability to high S_I values that are largely unrealistic. This kind of information was impossible to obtain so clearly by the standard Fisherian approaches. The rationale of choosing for S_I a prior more sophisticated than that employed for S_G, p₂, and G₀ was to determine with 95% confidence an interval as narrow as possible where S_I lies. The role played by this a priori knowledge is not hidden but transparent in terms of its effect on the results. For example, a strategy to obtain nonnegative confidence intervals by NLS could consist in reparameterizing the model with the aim of estimating the logarithm of the parameters. However, this would be obtained by approximating the likelihood as a function of the logarithm of the parameters, with a Gaussian distribution. This is questionable, because it would not be possible to evaluate the quality of such approximation.

As already mentioned, switching from NLS to Bayes estimation is not without a price, because the MCMC results are computationally more demanding to obtain. However, this strategy offers a new tool that could also become important in classification studies. In fact, it would be interesting to reanalyze large databases of epidemiological studies in which the standard NLS approach has produced S_I histograms with a likely artificial bimodal pattern (with one peak at S_I=0 and another peak at a positive S_I value), like that qualitatively displayed in Fig. 6, which are difficult to interpret. The value S_I=0, which Fisherian techniques return as the joint mode of a suitable objective function, can be a mathematical artifact, because it can be far from the minimum variance estimate provided by MCMC. Therefore, we believe that if the new strategy proposed in this paper were applied in epidemiological studies like those mentioned above, one would obtain a more reliable distribution of S_I estimates, e.g., changing the qualitative profile of Fig. 6A into that shown in Fig. 6B. In this way, more reliable information on NIDDM would be made available. Moreover, a reanalysis of databases such as GENNID (20), FUSION (26), and IRAS (16, 17), where the numbers of subjects are in the range of thousands, could also allow a refinement of our MCMC method, because a more accurate description of the a priori probability density function of S_I will be possible. It is clear, in fact, that the refined prior used in Fig. 5 provides better results than those obtained with the positivity prior and shown in Fig. 4. However, the refined prior still has a certain degree of roughness. A future evolution of the method could also envisage a population approach with the aim of estimating, jointly with model parameters, the common prior that has generated them.

View larger version (11K): [in this window] [in a new window]   Fig. 6.   A: a qualitative distribution of SI estimates that would be obtained in a study in which the SI=0 problem is present. On the x-axis, SI is measured in 104 min1/µU ml1. The y-axis shows the no. of SI estimates falling in equally spaced containers. B: qualitative distribution of SI estimates that would be obtained by the Markov chain Monte Carlo (MCMC) method in the same study.

Finally, it is worth stressing that important questions concerning possible problems introduced by the structure of the considered model still remain open also in the light of our MCMC results. In fact, one can clearly note that the a posteriori probability density functions of S_I of subjects 3 and 7 are quite close to the prior. As a consequence, the point estimates are close to the mean of the prior, equal to 1.67 10⁴ min¹/µU ml¹, suggesting that sometimes the model extracts little information on S_I from IVGTT data (in fact the posterior distribution approaches the prior as the data contain decreasing information about the parameters). Moreover, the presence of long tails in the a posteriori probability density of S_I of some subjects is a critical issue that will have to be clarified in future studies. Our Bayesian analysis, including reasonable a priori information on S_I, improves the estimation process, reducing this kind of problem. However, it is clear that no single parameter estimation technique, however sophisticated, can be a panacea, able to compensate possible structural defects of a model. So, further investigations are necessary to clarify whether the cause of these problems is the fact that the model is structurally inadequate for the estimation of a low S_I. We are inclined to believe that, in many of such critical situations, MCMC is really a more robust estimator than NLS or ML and can significantly improve the estimation process. However, we also believe that in the future it will also be worth exploring a different description of the same process, like the two-compartment model described in Refs. 7 and 8, which may overcome these problems. Maybe, after including some suitable a priori information, which is necessary because of a priori identifiability issues, the problems mentioned here may be solved by resorting to a maximum a posteriori estimator. In any case, it will be interesting to identify the two-compartment model by resorting to the same MCMC approach described here. In fact, even if more computationally demanding, this technique is likely to provide a clearer picture of the advantages of the two- vs. the one-compartment description.