Bayesian Identification of a Population Compartmental Model of C-Peptide Kinetics

Annals of Biomedical Engineering, Vol. 28, pp. 812–823, 2000 Printed in the USA. All rights reserved.

0090-6964/2000/28共7兲/812/12/$15.00 Copyright © 2000 Biomedical Engineering Society

Bayesian Identification of a Population Compartmental Model of C-Peptide Kinetics PAOLO MAGNI,1 RICCARDO BELLAZZI,1 GIOVANNI SPARACINO,2 and CLAUDIO COBELLI2 1

Dipartimento di Informatica e Sistemistica, Universita` degli Studi di Pavia via Ferrata 1, I-27100 Pavia, Italy and Dipartimento di Elettronica e Informatica, Universita` degli Studi di Padova, via Gradenigo 6a, I-35131 Padova, Italy

2

(Received 26 July 1999; accepted 3 July 2000)

Abstract—When models are used to measure or predict physiological variables and parameters in a given individual, the experiments needed are often complex and costly. A valuable solution for improving their cost effectiveness is represented by population models. A widely used population model in insulin secretion studies is the one proposed by Van Cauter et al. 共Diabetes 41:368–377, 1992兲, which determines the parameters of the two compartment model of C-peptide kinetics in a given individual from the knowledge of his/her age, sex, body surface area, and health condition 共i.e., normal, obese, diabetic兲. This population model was identified from the data of a large training set 共more than 200 subjects兲 via a deterministic approach. This approach, while sound in terms of providing a point estimate of C-peptide kinetic parameters in a given individual, does not provide a measure of their precision. In this paper, by employing the same training set of Van Cauter et al., we show that the identification of the population model into a Bayesian framework 共by using Markov chain Monte Carlo兲 allows, at the individual level, the estimation of point values of the C-peptide kinetic parameters together with their precision. A successful application of the methodology is illustrated in the estimation of C-peptide kinetic parameters of seven subjects 共not belonging to the training set used for the identification of the population model兲 for which reference values were available thanks to an independent identification experiment. © 2000 Biomedical Engineering Society. 关S0090-6964共00兲00907-3兴

partment model is usually adopted to describe CP kinetics.7,8 The model assumes that CP enters the system from the accessible compartment, where it is cleared from the system, and from which it distributes into a peripheral one 共Fig. 1兲. In order to provide the parameters of the model in a single individual, Polonsky and colleagues proposed in the 1980s21 a procedure entailing the execution of an input–output experiment: after having suppressed the spontaneous pancreatic secretion by means of a somatostatin infusion, an intravenous bolus of 共biosynthetic兲 CP is administered and plasma concentration samples are frequently collected 共for the sake of simplicity, this procedure will be hereafter referred to as the bolus experiment兲. Then, the two compartment model is fitted to the data to obtain individual estimates of model parameters. Obviously, this procedure is costly and labor intensive. In order to minimize both economical and ethical costs, an approach based on a population study was proposed in the early 1990’s by Van Cauter and co-workers.25 This approach takes advantage from the relatively small variability of CP kinetics among subjects in various physiopathological conditions and provides estimates of the individual parameters of the two compartment model from the knowledge of health status, sex, age and body surface area 共BSA兲 共hereafter referred to as the anthropometric parameters兲. Briefly, four deterministic linear regression models 共hereafter referred to as the population model兲 of short half-life (ts), long half-life (tl), amplitude fraction 共F兲 and volume of distribution (V) 共hereafter referred to as the noncompartmental parameters兲 were independently determined exploiting their values 共obtained from a bolus experiment兲 in a population of 111 normal subjects, 53 obese subjects, and 36 patients with noninsulin dependent diabetes mellitus with known anthropometric parameters. By using the population model, having obtained ts, tl, F, and V from age, sex, BSA, and health status of a given individual, the parameters of the two compartment model, i.e., k 01 , k 12 , k 21 and V 共hereafter referred to as the compartmental parameters兲 can be easily obtained by simple algebra. This

Keywords—C-peptide, Population model, Compartmental model, Bayes estimation, Markov chain Monte Carlo, Insulin, System identification

INTRODUCTION C-peptide 共CP兲 is cosecreted with insulin on an equimolar basis, exhibits linear kinetics in a large range of concentrations and, differently from insulin, is not extracted by the liver.9,21 Therefore, CP plays a key role in quantitative studies of the insulin system, e.g., plasma CP concentration data are used to estimate beta cell insulin secretion rate by deconvolution21,23 or to assess beta cell secretory indexes by the minimal model.24 In both cases the knowledge of the individual model of CP kinetics is required. In the literature, a linear two comAddress correspondence to Claudio Cobelli, Dipartimento di Elettronica e Informatica, Universita` degli Studi di Padova, via Gradenigo 6a, I-35131 Padova, Italy. Electronic mail: cobelli@[email protected]

812

Bayesian Identification of a Population Model of CP Kinetics

813

jects not belonging to the training set, for which reference estimates were independently obtained by a bolus experiment.22

PROBLEM FORMULATION FIGURE 1. The compartmental model of C-peptide kinetics.

approach is extensively adopted in both deconvolution3,5,15,16,20 and modeling24 studies on insulin secretion and its control by glucose. The approach adopted in Ref. 25 can be traced back to the class of two stage population modeling methods. In particular, in the first stage the individual kinetic model has been determined by nonlinear least squares and in the second stage the population model has been identified within a deterministic framework. This approach does not handle the existing sources of uncertainty. Therefore, a measure of the precision of the point estimates of the compartmental parameters k 01 , k 12 , k 21 and V in a given individual cannot be provided. However, the availability of this precision is crucial when one must make physiological inferences. For instance, a reliable confidence interval of the time course of insulin secretion reconstructed by deconvolution in an individual must take into account not only the uncertainty of CP plasma concentration data, but also that of the CP kinetic parameters used in the impulse response description. In order to arrive at a measure of the precision of the compartmental parameters k 01 , k 12 , k 21 and V in a given individual, it is necessary to embed the population model employed to determine noncompartmental parameters ts, tl, F and V into a probabilistic framework. This paper addresses within a Bayesian context the second step of the above described two stage procedure. In particular four new increasing-in-sophistication strategies are developed to identify the population model of Ref. 25. The Bayesian approach is, in general, more computer demanding than the deterministic one since stochastic simulation strategies, as Markov chain Monte Carlo methods, must be employed, but it is richer since it provides not only a point estimate of the model parameters but also their joint probability distribution. Our results on the identification of the population model show that for determining the joint distribution of the noncompartmental parameters it is of crucial importance to account for the uncertainty of the population model. After having determined the best strategy to identify the population model among the four possible choices, we evaluate its performance in the estimation of the individual compartmental parameters of seven normal sub-

The four CP noncompartmental parameters ts, tl, F, and V can be derived in each individual, without having to perform the bolus experiment, from the following four independent linear regressions:25 ts⫽m tsn ,

if subject is normal,

ts⫽m tso ,

if subject is obese,

ts⫽m tsd ,

if subject is diabetic,

F⫽m Fn ,

if subject is normal,

F⫽m Fo ,

if subject is obese,

F⫽m Fd ,

if subject is diabetic,

冧 冧

共1兲

共2兲

tl⫽a tl ⫹b tl * A, V⫽a Vm ⫹b Vm * BSA, V⫽a V f ⫹b V f * BSA,

if subject is a male, if subject is a female,

共3兲

冎

共4兲

where A is the age, expressed in years, BSA is expressed in m2 共BSA兲 is computed as: weight0.425 0.725 *height * 0.007 184, where the weight is expressed in kg and the height is expressed in cm兲 and m tsn , m tso , m tsd , m Fn , m Fo , m Fd , a tl , b tl , a Vm , b Vm , a V f , b V f are parameters of the population model. 共A subject is considered obese if his/her body weight is ⬎15% above ideal body weight.25兲 Let ␾ ts i , ␾ F i , ␾ tl i , ␾ V i be ts, F, tl, V in the ith subject, and ␪ ts ⫽ 关 m tsn ,m tso ,m tsd 兴 T , ␪ F T ⫽ 关 m Fn ,m Fo ,m Fd 兴 , ␪ tl ⫽ 关 a tl ,b tl 兴 T , and ␪V ⫽ 关 a Vm ,b Vm ,a V f ,b V f 兴 T be the parameter vectors of the population model. Then,

␾ ts i ⫽U ts i ␪ ts ⫹ ␩ ts i ,

共5兲

␾ F i ⫽U F i ␪ F ⫹ ␩ F i ,

共6兲

␾ tl i ⫽U tl i ␪ tl ⫹ ␩ tl i ,

共7兲

␾ V i ⫽U V i ␪ V ⫹ ␩ V i ,

共8兲

814

MAGNI et al.

where ␩ ts i , ␩ F i , ␩ tl i , ␩ V i , are the differences between the available data 共i.e., ␾兲 and the values predicted by the model 共i.e., U ␪ ); U ts i , U F i , U tl i , U V i are suitable matrices containing age, sex, BSA, and health status of the ith subject. It is easy to see that

U ts i ⫽

U Fi⫽

冋再 再 冋再 再

再 册 再 册

1

if normal 1

if obese

0

otherwise, 0

otherwise,

1

if diabetic

0

otherwise

,

1

if normal 1

if obese

0

otherwise, 0

otherwise,

1 0

冋再 再

再 再

BSA

k 21⫽

1

if male

0

otherwise, 0

otherwise,

1

if female BSA

if female

0

otherwise, 0

otherwise

ln共 2 兲 ln共 2 兲 1 , ts tl k 12

共12兲

ln共 2 兲 ln共 2 兲 ⫹ ⫺k 12⫺k 01 . ts tl

共13兲

F 1⫺F ⫹ , tl ts

册

IDENTIFICATION OF THE POPULATION MODEL .

␾ i ⫽U i ␪ ⫹ ␩ i . Assuming that both noncompartmental and anthropometric parameters are available in N subjects 共the training set兲, the whole model for each of the four regressions can be written as

␾ ⫽U ␪ ⫹ ␩ ,

共9兲

where ␾ ⫽ 关 ␾ 1 , ␾ 2 ,..., ␾ N 兴 , U⫽ 关 U 1 ,U 2 ,...,U N 兴 , and ␩ ⫽ 关 ␩ 1 , ␩ 2 ,..., ␩ N 兴 T . For the sake of reasoning it is useful to split ␩ in two components, assumed to be additive

␩ ⫽ ␩ 1⫹ ␩ 2 ,

冉

if male

For the sake of simplicity, in the rest of the paper we will adopt for the generic regression model among Eqs. 共5兲–共8兲 the following notation:

T

共11兲

k 01⫽

,

U tl i ⫽ 关 1, A 兴

U Vi⫽

冊

k 12⫽ln共 2 兲

if diabetic otherwise

共the so-called random effects兲 and will be referred to as ‘‘modeling error.’’ Once the identification of the population model is performed, i.e., the vector ␪ is derived for each regression model as discussed in the next section, one can move to the individual level, where the noncompartmental parameters of a generic subject not included in the training set, denoted by ␾ ⫹ , can be obtained from his/ her anthropometric parameters U ⫹ . Finally, the compartmental parameters can be derived by the following algebraic equations:

T

共10兲

where ␩ 1 accounts for the uncertainty of the data and will be referred to as ‘‘data error.’’ ␩ 2 accounts for the fact that subjects with the same anthropometric parameters can have different noncompartmental parameters

The training-set available for the identification of the CP population model is a slightly updated version of that used in Ref. 25. It consists of 110 normal subjects 共71 male, 39 female兲, 52 obese subjects 共18 male, 34 female兲, and 45 diabetic patients 共25 male, 20 female兲. For each subject both the anthropometric parameters and the estimates of the noncompartmental parameters are available 共courtesy of Dr. Polonsky兲. The estimates of noncompartmental parameters were obtained in Ref. 25 by fitting the model of CP kinetics against experimental data obtained, for each subject, by measuring the CP plasma concentration after a bolus administration.

Model 0: Ignoring All the Sources of Error When no sources of uncertainty are considered in the population model, the estimation problem Eq. 共9兲 can be solved in a deterministic framework by using ordinary least squares. This was the approach adopted in Ref. 25 and hereafter will thus be referred to as Model 0 共M0兲. M0 parameter estimates are

␪ˆ ⫽ 共 U T U 兲 ⫺1 U T ␾ .

共14兲

Having obtained ␪ˆ by applying Eq. 共14兲 to each regression model Eqs. 共5兲–共8兲, the individual noncompartmental parameters of a new subject can be calculated as

Bayesian Identification of a Population Model of CP Kinetics

␾ ⫹ ⫽U ⫹ ␪ˆ .

共15兲

In this context, the point estimate of the population model parameters is derived with no information on their precision. As a consequence, the determination of the individual CP noncompartmental parameters through Eq. 共15兲 共and thus that of the compartmental parameters兲 cannot take into account the interindividual variability and a measure of the precision of the individual estimates is not available. In order to derive the precision together with the point estimate of both the population and the individual parameters, the problem can be conveniently formulated into a Bayesian framework, where all model variables 共i.e., data, regression parameters, errors兲 are stochastic variables described by probability distributions. In the following, four population models based on different statistical hypotheses will be considered. The first three models 共M1, M2, M3兲 consider the four regressions Eqs. 共5兲–共8兲 independently 关the major difference among these models lies in the probabilistic assumptions on ␩ in Eq. 共9兲兴, so that the estimation problem Eq. 共9兲 can be solved by considering four distinct estimation problems. In contrast, the fourth model 共M4兲 is able to take into account the correlations between the four regression models, so that all the regression parameters have to be jointly estimated. Model 1: Data Error Only In order to account for ␩, the simplest situation is that which assumes no interindividual variability. In other words, the only source of uncertainty kept into account by M1 is that of the estimates of the noncompartmental parameters of the training set. ␩ i are assumed independent from each other, and in particular

␩ ⬃N 共 0,⌿ 兲 ,

共16兲

where N(•,•) is Normal distribution and ⌿ is a known matrix containing the variance of the data ␾. 关Matrix ⌿ quantifies the uncertainty of the noncompartmental parameters ␾ obtained from the bolus experiment data 共affected by a 4% measurement error CV兲 and is here determined as the inverse of the Fisher information matrix associated to the bolus experiment.兴 Moreover, M1 assumes that no information is available on ␪. In symbols

␪ ⬃Uni f 关 ⫺⬁,⬁ 兴 ,

共17兲

where Unif is the uniform distribution. 关An alternative choice for the prior distribution of ␪ could be the Normal distribution with a large variance in order to obtain a

815

poor informative distribution. However, when the variance of the prior distribution tends to infinity, the posterior distribution tends to Eq. 共18兲.兴 Given the stochastic model Eqs. 共16兲–共17兲, the posterior distribution of ␪ is calculated as p 共 ␪ 兩 ␾ 兲 ⬀p 共 ␾ 兩 ␪ 兲 p 共 ␪ 兲 ⫽N 共 ␮ ,⍀ 兲

共18兲

with

␮ ⫽ 共 U T ⌿ ⫺1 U 兲 ⫺1 U T ⌿ ⫺1 ␾ , ⍀⫽ 共 U T ⌿ ⫺1 U 兲 ⫺1 . It is interesting to note that the Bayesian point estimate 共i.e., the mean of the posterior distribution兲 of this model is the same as that of weighted least squares. Then, the probabilistic model to predict individual noncompartmental parameters, given those of the training set, becomes T p 共 ␾ ⫹ 兩 ␾ 兲 ⫽N 共 U ⫹ ␮ ,U ⫹ ⍀U ⫹ 兲

共19兲

being ␾ ⫹ ⫽U ⫹ ␪ . This approach, differently from the previous one 共M0兲, is able to consider in the calculation of the individual parameters Eq. 共19兲 the data error ␩ 1 . However, in population studies it is of crucial importance to assess the role of the random effects, here the modeling error ␩2 .

Model 2: Data and Modeling Errors In this model, called M2, we consider the fact that ␩ in Eq. 共9兲 is due to both data and modeling errors. Having assumed in Eq. 共10兲 the statistical independence of ␩ 1 and ␩ 2 , one has

␩ ⬃N 共 0,⌿⫹ ␴ 2 I 兲 ,

共20兲

where ⌿ is the matrix introduced in Eq. 共16兲 and ␴ 2 is an unknown parameter to be estimated. 共For generality, in each regression, we consider a different variance parameter ␴ 2 that can be conveniently denoted with 2 ␴ ts , ␴ F2 , ␴ 2tl or ␴ V2 .) In order to completely describe this stochastic model, we have to make a choice on prior distributions. In particular we assume that

␪ ⬃N 共 ␪ 0 ,⌺ 0 兲 ,

共21兲

816

MAGNI et al.

␴ ⫺2 ⬃⌫ 共 ␥ 1 , ␥ 2 兲 ,

共22兲

where ␪ 0 , ⌺ 0 , ␥ 1 , ␥ 2 are fixed parameters. ⌫( ␥ 1 , ␥ 2 ) is the gamma distribution with mean ␥ 1 / ␥ 2 and variance ␥ 1 / ␥ 22 . In this case the joint posterior distribution of the stochastic parameters 兵 ␪ , ␴ 2 其 cannot be derived in a closed form, and it is necessary to resort to a stochastic simulation strategy known as Markov chain Monte Carlo 共MCMC兲. MCMC methods are based on two steps: Markov chain generation and Monte Carlo integration.13 Briefly, by sampling from suitable probability distributions, it is possible to generate a Markov chain that converges 共in distribution兲 to the target distribution, i.e., the distribution to be derived. MCMC methods differ from each other in the way the Markov chain is created, even if all the different strategies proposed in the literature can be traced back to the Metropolis–Hastings algorithm.14,19 For the sake of readability, we report in the Appendix the technical details on the specific sampling scheme we developed to cope with this model identification problem. Having obtained the sample posterior distribution of the population parameters p( ␪ , ␴ 2 兩 ␾ ), it is possible to derive the distribution of the noncompartmental parameters of the new subject by computing p 共 ␾ ⫹兩 ␾ 兲 ⫽

冕

␪,␴2

p 共 ␾ ⫹兩 ␪ , ␴ 2, ␾ 兲 p 共 ␪ , ␴ 2兩 ␾ 兲 d ␪ d ␴ 2, 共23兲

where

␩ 2 is much greater than the data error ␩ 1 共in our hands about 50 times兲, so that the latter can be neglected. These remarks suggest reformulating M2 as described in the next section. Model 3: Modeling Error Only If we neglect the data error ␩ 1 we obtain model 3 共M3兲, where

␩ ⬃N 共 0,␴ 2 I 兲

with ␴ 2 being an unknown parameter to be estimated from the data. One can use the following assumptions on the probability distributions of the population parameters:

␪ ⬃U 关 ⫺⬁,⬁ 兴 ,

共25兲

␴ ⫺2 ⬃⌫ 共 ␥ 1 , ␥ 2 兲 ,

共26兲

where ␥ 1 , ␥ 2 are fixed parameters. In this case the posterior distribution of the stochastic parameters 兵 ␪ , ␴ 2 其 can be derived in closed form. In fact, via some computations p 共 ␪ , ␴ 2 兩 ␾ 兲 ⬀p 共 ␾ 兩 ␪ , ␴ 2 兲 p 共 ␪ 兲 p 共 ␴ 2 兲 ⬀⌫ ␴ ⫺2 共 a,b 兲 N ␪ 共 ␮ ,⍀ 兲 共27兲 with a⫽ ␥ 1 ⫹

p 共 ␾ ⫹ 兩 ␪ , ␴ 2 , ␾ 兲 ⫽N 共 U ⫹ ␪ , ␴ 2 兲 . Obviously, since p( ␪ , ␴ 2 兩 ␾ ) is in sample form, p( ␾ ⫹ 兩 ␾ ) must be obtained in sample form too by using, e.g., the following forward sampling procedure: 共1兲 randomly extract one of the available samples of p( ␪ , ␴ 2 兩 ␾ ); 共2兲 extract one sample from p( ␾ ⫹ 兩 ␪ , ␴ 2 , ␾ ); and 共3兲 repeat steps 共1兲 and 共2兲 until the calculation of expectations becomes reliable. In conclusion, it is important to make two remarks on this model. First, the computational burden required to draw samples from p( ␪ , ␴ 2 兩 ␾ ) is significant. However, this step has to be performed only once, during the population model identification. In contrast, the Monte Carlo simulation required to obtain p( ␾ ⫹ 兩 ␾ ) for each new subject is not time expensive. Second, in the CP population model identification problem, as it can be seen from Figs. 2 and 3 of Ref. 25, the modeling error

共24兲

b⫽ ␥ 2 ⫹

共 N⫺q⫺2 兲 , 2

共 y⫺U ␮ 兲 T 共 y⫺U ␮ 兲 , 2

␮ ⫽ 共 U T U 兲 ⫺1 U T y, ⍀⫽ ␴ 2 共 U T U 兲 ⫺1 , where q is the number of the regression parameters 共length of ␪兲, ⌫ ␴ ⫺2 (•,•) is a gamma distribution of the variable ␴ ⫺2 , and N ␪ (•,•) is a normal distribution of the variable ␪. The distribution p( ␾ ⫹ 兩 ␾ ) of the noncompartmental parameters in a new subject can still be obtained by Eq. 共23兲. Although the analytic expression of the posterior distribution p( ␪ , ␴ 2 兩 ␾ ) is available, only the sample distribution of p( ␾ ⫹ 兩 ␾ ) can be computed, by using for example the following forward sampling scheme: 共1兲 extract one sample for ␴ 2 from ⌫ ␴ ⫺2 (a,b);

Bayesian Identification of a Population Model of CP Kinetics

共2兲 extract one sample for ␪ from N ␪ ( ␮ ,⍀); 共3兲 extract one sample for ␾ ⫹ from p( ␾ ⫹ 兩 ␪ , ␴ 2 , ␾ ); and 共4兲 repeat steps 共1兲, 共2兲 and 共3兲. In M3, the slight approximation introduced by neglecting the data error in ␩ has allowed us to greatly simplify the computational burden involved in the identification of the population model. Moreover, it is not required to store samples of the posterior distribution of population parameters p( ␪ , ␴ 2 兩 ␾ ), but only the value of the parameters a, b, ␮ and ⍀. Therefore, this model appears as the most convenient to tackle our problem when the four regression models are treated separately. However, since the noncompartmental parameters are likely to be correlated, it is worthwhile to investigate the performance of a model where correlation among the modeling errors is explicitly taken into account.

are assumed to be independent from each other, so that we have ˜ 兲, ˜␩ ⬃N 共 0,⌺ with

冋

˜ ⫽U ˜ ˜␪ ⫹ ˜␩ , ␾ i i i

共28兲

˜ ⫽ 关 ␾ , ␾ , ␾ , ␾ 兴 T is the vector of the nonwhere ␾ i ts i Fi tl i Vi compartmental parameters for the ith subject, ˜␪ ⫽ 关 ␪ ts , ␪ F , ␪ tl , ␪ V 兴 T is the vector of population parameters, ˜␩ i ⫽ 关 ␩ ts i , ␩ F i , ␩ tl i , ␩ V i 兴 T is the vector of errors and ˜ is a 共4⫻12兲 matrix containing age, sex, BSA and U i health status of the ith subject. It is easy to see that

˜ ⫽ U i

冋

U ts i

0

0

0

0

U Fi

0

0

0

0

U tl i

0

0

0

0

U Vi

册

¯ 

¯

0

共30兲

册

] , ⌺

˜ is a (4N⫻4N) matrix and ⌺ is a 共4⫻4兲 matrix where ⌺ to be estimated from the data. The distributions needed to fully specify the stochastic model are ˜␪ ⬃N 共 ˜␪ ,⌺ ˜ 兲, 0 0

共31兲

˜ 兲, ⌺ ⫺1 ⬃W 共 ˜␳ ,R

共32兲

˜ , ˜␳ , ˜R fixed parameters and W is the with ˜␪ 0 , ⌺ 0 11 Wishart distribution with mean ␳ R. Since also in this case the posterior distribution of the stochastic parameters 兵˜␪ ,⌺ 其 cannot be derived in closed form, we again resort to an MCMC method 共details in the Appendix兲. Having obtained the sample posterior distribution of ˜ ), it is possible to the population parameters p(˜␪ ,⌺ 兩 ␾ derive the distribution of the noncompartmental individual parameters: ˜ 兩␾ ˜ 兲⫽ p共 ␾ ⫹

冕

˜␪,⌺

˜ 兩˜␪ ,⌺, ␾ ˜ 兲 p 共 ˜␪ ,⌺ 兩 ␾ ˜ 兲 d˜␪ d⌺, p共 ␾ ⫹ 共33兲

where ˜ 兩˜␪ ,⌺, ␾ ˜ 兲 ⫽N 共 U ˜ ˜␪ ,⌺ 兲 . p共 ␾ ⫹ ⫹ ˜ ), it is possible From the sample distribution p(˜␪ ,⌺ 兩 ␾ ˜ 兩␾ ˜ ), by exto derive the sample distribution of p( ␾ ⫹ ploiting the following forward sampling strategy:

.

Adopting a rationale similar to that employed to arrive at Eq. 共9兲, the population model can be written as ˜ ⫽U ˜ ˜␪ ⫹ ˜␩ , ␾

⌺

˜⫽ ] ⌺ 0

Model 4: Correlation Among Modeling Errors Model 4 共M4兲 is a generalization of M3, which takes into account correlations among modeling errors of the four regressions Eqs. 共5兲–共8兲 in the same subject. To this aim, it is impossible to consider separately the four estimation problems, and it is necessary to group the regressions in a single model, which is then used for a joint estimate of all regression parameters 兵 ␪ ts , ␪ F , ␪ tl , ␪ V 其 and error variances. It is possible to write, for the ith subject

817

共29兲

˜ ⫽关␾ ˜ ,␾ ˜ ,..., ␾ ˜ 兴 T, U ˜ ⫽关U ˜ ,U ˜ ,...,U ˜ 兴 T , and where ␾ 1 2 N 1 2 N ˜␩ ⫽ 关 ˜␩ 1 , ˜␩ 2 ,..., ˜␩ N 兴 T . Errors related to different subjects

共1兲 randomly extract one of the available samples of ˜ ); p(˜␪ ,⌺ 兩 ␾ ˜ 兩˜␪ ,⌺, ␾ ˜ ); and 共2兲 extract one sample from p( ␾ ⫹ 共3兲 repeat steps 共1兲 and 共2兲. RESULTS In order to completely specify M2, M3, and M4, it is necessary to assign values to the parameters involved in the a priori probabilistic model. In particular:

818

MAGNI et al.

TABLE 1. Point estimate of the population model parameters obtained by applying both deterministic „M0… and Bayesian approaches „M1, M2, M3, M4… on the training data set.

␪

M0

M1

M2

M3

M4

m tsn m tso m tsd m Fn m Fo m Fd a tl b tl a Vm b Vm a Vf b Vf

5.00 4.55 4.59 0.764 0.782 0.780 27.79 0.177 0.491 1.99 1.55 1.42

4.67 4.12 4.16 0.781 0.794 0.797 25.56 0.098 0.587 1.88 1.44 1.37

4.99 4.54 4.58 0.765 0.782 0.780 27.80 0.176 0.795 1.46 1.84 1.46

5.00 4.55 4.60 0.764 0.782 0.779 27.81 0.177 0.477 1.55 2.00 1.41

4.99 4.50 4.70 0.767 0.781 0.778 26.70 0.209 0.417 0.89 2.04 1.77

TABLE 2. Point estimate of the population model error variance, obtained by applying Bayesian approaches M2, M3, and M4 on the training data set. Modeling error variance M2 M3 ‘ M4

␴ 2ts ⫽1.20, ␴ 2F ⫽0.0015, ␴ tl2 ⫽32.82, ␴ V2 ⫽0.678 ␴ 2ts ⫽1.35, ␴ 2F ⫽0.0017, ␴ tl2 ⫽34.33, ␴ V2 ⫽0.740 1.30 0.0062 3.24 0.59

冋

⌺⫽

Corr共 ⌺ 兲 ⫽ 共i兲

共ii兲 共iii兲

For M2, in Eqs. 共21兲 and 共22兲, we have assumed ␪ 0 ⫽ 关 5 5 5 兴 T for the first regression, ␪ 0 ⫽ 关 1 1 1 兴 T for the second regression, ␪ 0 ⫽ 关 30 1 兴 T for the third regression, and ␪ 0 ⫽ 关 1 1 1 1 兴 T for the last one; ⌺ 0 was taken as a diagonal matrix with ␪ 0 square elements on the diagonal, ␥ 1 ⫽0.1, ␥ 2 ⫽2e7. For M3 in Eq. 共26兲, the values of ␥ 1 , ␥ 2 specified for M2 have been used. For M4 in Eqs. 共31兲 and 共32兲, we have chosen ˜␪ ⫽ 关 5 5 5 1 1 1 30 1 1 1 1 1 兴 T , ⌺ ˜ as a diagonal 0 0 matrix with ˜␪ 0 square elements on the diagonal, ␳⫽10 and ˜R ⫽ ␳ ⫺1 (0.01* diag(关5 1 30 4 兴 )) ⫺1 , where diag is the diagonal matrix operator.

These parameters allow us to obtain prior distributions centered around reasonable values 共derived from Ref. 25兲, even if ‘‘sufficiently flat.’’ In Table 1 we report the point estimates of the population parameters. Although, in general, the point estimates do not change significantly with the strategy adopted, one must keep in mind that the Bayesian approaches provide much richer information given by the posterior distribution of all parameters. In Table 2 we report the point estimates of the variance of the modeling error Eq. 共20兲 for M2, Eq. 共24兲 for M3, and Eq. 共30兲 for M4. M3 has modeling error variances slightly larger than those of M2. This is not surprising, since in M2 also a data error component is present 关see Eq. 共20兲兴. The diagonal elements of the variance matrix ⌺ of M4 are comparable to the variance values obtained for M2 and M3. The elements outside the diagonal are more easily interpreted by computing the correspondent correlation matrix:

冋

0.0022

0.071

⫺0.0055

3.24

0.071

32.94

1.89

0.59

⫺0.0055

1.89

0.71

0.0062

1

0.12

0.12

1

0.50

0.27

0.61 ⫺0.14

0.50

0.61

0.27 ⫺0.14 1

0.39

0.39

1

册

.

册

共34兲

This matrix highlights a significant correlation between model variables thus making the assumption of independence between regression models not tenable. For instance, correlation is found between the regression errors of V and ts 共0.61兲, between ts and tl 共0.50兲, or between V and tl 共0.39兲. In order to select the best model among M1–M4, we have adopted the Schwarz–Bayesian information criterion, which gives the following scores: ⫺8.09e4 for M1, ⫺1.85e3 for M2, ⫺1.86e3 for M3, and ⫺1.7e3 for M4. 关The best model is the one having the highest score computed using the formula BIC⫽2l⫺p log n, where l is the expected value of the log likelihood, p is the number of parameters in the model, and n is the number of data points. In our case n⫽207* 4 and p is 12 共sum of size of ␪ ts , ␪ F , ␪ tl , ␪ V ), 16 共12 parameters for ␪ s plus 2 ␴ ts , ␴ F2 , ␴ 2tl , ␴ V2 ), 16 共as before兲, 22 共12 parameters for ˜␪ plus 10 parameters for ⌺兲 for M1, M2, M3, and M4, respectively.兴 The superiority of M4 is not surprising since M4 is the only model that relaxes the assumption on the incorrelation of the four regression models, the criticality of which is suggested by the nondiagonal elements of the matrix Corr共⌺兲. INDIVIDUAL PARAMETERS PREDICTED BY THE POPULATION MODEL: EVALUATION AGAINST ‘‘TRUE’’ VALUES The ultimate goal of the CP population model is to provide the values of the compartmental parameters in a given subject from the knowledge of his/her anthropometric parameters. Therefore, it is helpful to compare the values of the compartmental parameters obtained by the population approach with those provided through an independent experiment, in which the compartmental

Bayesian Identification of a Population Model of CP Kinetics TABLE 3. Anthropometric parameters of the 7 normal subjects used to evaluate the prediction of the stochastic population model. Subject No.

Sex

Age (yr)

Height (cm)

Weight (kg)

1 2 3 4 5 6 7

Male Male Female Female Female Male Male

25 22 33 19 28 21 23

181.8 185 158 168.8 164 185.4 165

70.7 83.6 52 67 56.4 94.3 68.2

model parameters have been obtained by identifying the model from samples following a bolus administration. Here we consider a group of seven subjects not belonging to the training set.22 For this set of subjects, differently from the ones in the training set, the ‘‘raw’’ experimental data coming from the bolus experiment are available, as well as the anthropometric parameters reported in Table 3. Therefore, we compare the compartmental parameters obtained from the bolus experiment with those predicted by using M0, which is the model presently used in the literature, and M4. 关The joint prob˜ ) is obtained in ability distribution p(k 12 ,k 01 ,k 21 ,V 兩 ␾ M4 solving the algebraic equation Eqs. 共11兲–共13兲 for ˜ 兩␾ ˜ ), derived each sample of the distribution p( ␾ ⫹ through Eq. 共33兲.兴 Figure 2 shows for subject No. 1 the point estimates of the compartmental parameters together with their 95% confidence intervals provided from the bolus experiment; the 共point兲 estimates obtained by M0; and the point es-

819

timates and the confidence intervals obtained by M4. In three over four cases 共i.e., all except k 21) the M0 point estimates lie outside the 95% confidence intervals of the bolus estimates. In other words, the M0 estimates have a large probability of being incorrect. In contrast, the new approach 共M4兲 provides an interval that contains the 共bolus兲 band to which the ‘‘true value’’ belongs. In Table 4 the same information reported in graphical form in Fig. 2 for subject No. 1 is given for all the subjects. Comments similar to those drawn above can be made. Results clearly demonstrate the importance of quantifying the uncertainty of the individual compartmental parameters obtained by the population approach. Obviously, inferences made by using these parameters must account for this uncertainty. To this aim, it is necessary to consider the joint posterior distribution ˜ ), which allows one to account for p(k 12 ,k 01 ,k 21 ,V 兩 ␾ the correlation among parameters 共this information is not represented in Table 4, in which the confidence intervals reflect only the marginal distribution of each parameter兲. All the algorithms were implemented in MATLAB 共The MathWorks, Inc., Natick, MA兲. In order to give an idea of computational times, about 10 min were required to identify M4 on a Sun Sparc Station Ultra from the data of the training set, while the estimation of the individual parameters in a new subject 共not belonging to the training set兲 by M4 only required a few seconds. CONCLUSIONS When models are used to measure/predict variables and parameters in a given individual, a crucial problem

FIGURE 2. CP compartmental parameters for subject No. 1. The point estimates are denoted by a line and the 95% confidence intervals by a gray zone.

820

MAGNI et al. TABLE 4. Estimates of the CP two compartment model parameters in the 7 subjects. The parameters obtained from the bolus experiment „Bolus… are compared with those obtained by M0 and M4 population approaches. The standard deviation is in parenthesis, whereas the 95% confidence interval is in square brackets. Subject No. Bolus 1

M0 M4 Bolus

2

M0 M4 Bolus

3

M0 M4 Bolus

4

M0 M4 Bolus

5

M0 M4 Bolus

6

M0 M4 Bolus

7

M0 M4

k 12* 102 (min⫺1)

k 01* 102 (min⫺1)

k 21* 102 (min⫺1)

V (l)

6.17 (0.66) [4.96, 7.65] 4.91 5.18 (1.41) [3.23, 8.68]

8.45 (0.72) [7.14, 10.05] 6.07 6.37 (1.28) [4.44, 9.34]

6.92 (1.43) [4.42, 10.03] 5.03 5.42 (2.43) [2.69, 11.31]

3.43 (0.33) [2.80, 4.11] 4.30 4.30 (0.86) [2.60, 6.00]

5.86 (0.71) [4.61, 7.45] 4.94 5.26 (1.52) [3.32, 8.88]

8.51 (0.82) [6.99, 10.02] 6.14 6.46 (1.28) [4.53, 9.47]

7.25 (1.62) [4.47, 10.82] 4.97 5.45 (2.64) [2.58, 11.78]

3.25 (0.35) [2.60, 4.01] 4.63 4.64 (0.84) [2.93, 6.28]

5.68 (0.65) [4.46, 7.06] 4.84 5.15 (1.53) [3.19, 8.69]

6.54 (0.57) [5.45, 7.70] 5.90 6.17 (1.21) [4.41, 8.95]

7.16 (1.51) [4.43, 10.45] 5.18 5.64 (2.59) [2.75, 12.06]

3.29 (0.31) [2.73, 3.95] 3.70 3.56 (0.86) [1.89, 5.21]

6.40 (0.63) [5.29, 7.75] 4.97 5.25 (1.40) [3.28, 8.79]

8.35 (0.75) [6.97, 9.92] 6.20 6.53 (1.31) [4.55, 9.72]

9.05 (1.72) [6.02, 12.8] 4.91 5.25 (2.28) [2.56, 11.24]

2.85 (0.29) [2.31, 3.45] 4.06 4.03 (0.85) [2.30, 5.73]

6.46 (0.74) [5.18, 8.12] 4.89 5.21 (1.48) [3.20, 8.84]

6.08 (0.53) [5.10, 7.18] 6.00 6.31 (1.25) [4.43, 9.23]

9.27 (1.88) [6.09, 13.44] 5.09 5.56 (2.54) [2.73, 11.51]

3.38 (0.32) [2.80, 4.04] 3.83 3.73 (0.84) [2.07, 5.37]

7.63 (0.93) [6.01, 9.68] 4.95 5.27 (1.47) [3.22, 8.79]

8.83 (0.87) [7.33, 10.67] 6.16 6.49 (1.27) [4.61, 9.46]

10.90 (2.39) [6.94, 16.19] 4.95 5.38 (2.51) [2.54, 11.68]

3.90 (0.43) [3.10, 4.77] 4.85 4.87 (0.84) [3.21, 6.50]

5.41 (0.38) [4.59, 6.08] 4.93 5.27 (1.51) [3.20, 8.74]

5.16 (0.20) [4.74, 5.52] 6.11 6.49 (1.33) [4.51, 9.68]

3.44 (0.37) [2.54, 3.93] 4.99 5.46 (2.63) [2.54, 11.52]

3.95 (0.17) [3.66, 4.33] 3.98 3.96 (0.85) [2.29, 5.62]

is related to the ethical, technological, and economical costs of the often complex experiments involved. Population models may be of great help in improving cost effectiveness, and this makes them appealing. For instance, a widely used population model allows the determination of the four parameters of a two compartment model of CP kinetics in a given subject from the knowledge of his/her age, sex, BSA, and health condition, without having to perform an ad hoc bolus experiment. When physiologic parameters in an individual are obtained through a population approach, it would also be highly desirable to have a measure of their precision. However, when population approaches are ‘‘determinis-

tic’’ it is impossible to arrive at the precision of the parameter estimates in the single individual. Various approaches to population modeling, which take into account stochastic components, are available in the literature, e.g., in pharmacokinetics, and most of them can be traced back to the so called global approaches.2,6,10,18,26,27 These approaches jointly consider the individual and the population levels and thus make the propagation of the uncertainty from data to population parameters possible. Unfortunately, a global approach is not suitable to our problem because the experimental data at the individual level 共i.e., the CP decay curves of Ref. 25兲 are not available. It is worth noting

Bayesian Identification of a Population Model of CP Kinetics

821

FIGURE 3. CP unit impulse response in subject No. 7. The solid line is the prediction when the bolus data are used. The dashed line is the estimate curve obtained by using both population model M3 „considering only modeling error… and M4 „considering correlation among modeling errors…. The light gray region is the 95% confidence interval obtained with M4, while the light gray plus dark gray region gives the 95% confidence interval obtained with M3.

however that situations similar to ours can arise frequently in physiological system studies, e.g., from the literature one is only able to build a data base of ‘‘aggregated’’ parameters by gathering information from studies published by different research groups.1 In this paper we have investigated four strategies to identify the widely used population model of the CP kinetics, proposed in Ref. 25 within a Bayesian framework taking into account the existing sources of uncertainty, i.e., data and modeling error 共we felt it was outside the scope of our work to explore possible new population models of the CP kinetics兲. Our results, obtained on a slightly updated version of the training set of Ref. 25, show that data error can be neglected in comparison with modeling error and that accounting for the correlation between the regression models of the noncompartmental parameters is of key importance. The stochastic formulation of the population model allows us to derive the compartmental parameters together with their precision in the single individual, as shown in our evaluation study on seven subjects 共not belonging to the training set兲 for which reference parameters were available from a bolus experiment. The importance of explicitly accounting for the correlation among the regression models of the noncompartmental parameters emerges even more clearly if we consider the following example. The CP noncompartmental parameters can be used to reconstruct the unit impulse response

the unit impulse response obtained from the bolus experiment22 and the point estimate with 95% confidence bands provided by M3 and M4. This is evident from M3, ignoring the correlation between the errors in the four regressions results in overly pessimistic confidence intervals. Similar considerations can be obtained for each of the other subjects. In conclusion, the paper shows the importance of obtaining from the population model not only the point estimates of the compartmental parameters in a given individual, but also their precision. This precision should be explicitly considered when physiological inferences are made using the individual kinetic parameters. For instance, a reliable confidence interval of insulin secretion rate reconstructed by deconvolution must take into account not only the uncertainty of the CP plasma concentration data, but also that of the individual kinetic parameters.

y 共 t 兲 ⫽1/V 共 Fe ⫺ln共 2 兲 /ts ⫹ 共 1⫺F 兲 e ⫺ln共 2 兲 /tl 兲 .

APPENDIX

Figure 3 shows, for a representative subject 共No. 7兲,

ACKNOWLEDGMENTS We would like to thank Dr. K. S. Polonsky 共Department of Medicine, University of Chicago, IL兲 for having made available to us the CP data of Ref. 22 and the data base of Ref. 25. This work was in part supported by NIH Grant Nos. RR-11095 and RR-12609, by a University of Padova grant 共STIM-PET兲 to G. S., and by a University of Pavia grant to P. M.

Here we discuss the technical details on the MCMC scheme built for the identification of the M2 and M4

MAGNI et al.

822

population models. The most popular MCMC scheme is the Gibbs sampling12 that allows us to extract samples for each variable from full conditional distributions, i.e., the probability distribution of a variable given all others. Similarly, our MCMC schemes require us to draw samples directly from the full conditionals when it is possible.

˜ 兲 ⬀p 共 ˜␪ , ␾ ˜ 兩 ⌺ 兲 ⫽ p 共 ˜␪ 兲 p 共 ␾ ˜ , 兩˜␪ ,⌺ 兲 ⫽N 共 ␮ ,⍀ 兲 , p 共 ˜␪ 兩 ⌺, ␾ 共A3兲 ˜ 兲 ⬀p 共 ⌺, ␾ ˜ 兩˜␪ 兲 ⫽ p 共 ⌺ 兲 p 共 ␾ ˜ , 兩˜␪ ,⌺ 兲 ⫽W 共 ␳ ,R 兲 , p 共 ⌺ 兩˜␪ , ␾ 1 1 共A4兲 with

冉

˜ ⫺1 ⫹ ␮⫽ ⌺ 兺 U˜ Ti ⌺ ⫺1 U˜ i 0

Model 2 In this model the full conditional distributions are

˜ ⫺1 ⫹ ⍀⫽ ⌺ 0

冉 冏 冊 冉 冏冊 冉 冊 1 1 ,␾ ␪ 2 ␪ , ␾ ⬀p ␴ ␴2 ⫽p

冉

1 p共 ␾兩␪,␴2兲 ␴2

R 1 ⫽ ˜R ⫺1 ⫹

˜ ˜ T ⌺ ⫺1 ␾ U i i

兺i U˜ Ti ⌺ ⫺1 U˜ i

冊

冊

⫺1

,

⫺1

,

⫻ 共 ⌿⫹ ␴ 2 I 兲 ⫺1 共 ␾ ⫺U T ␪ 兲兲 ,

共A2兲

where T 2 ⫺1 ␮ ⫽ 共 ⌺ ⫺1 U 兲 ⫺1 0 ⫹U 共 ⌿⫹ ␴ I 兲

⫻ 共 U T 共 ⌿⫹ ␴ 2 I 兲 ⫺1 ␾ ⫹⌺ ⫺1 ␪ ␪0兲, T 2 ⫺1 ⍀⫽ 共 ⌺ ⫺1 U 兲 ⫺1 . 0 ⫹U 共 ⌿⫹ ␴ I 兲

Unfortunately, in this case we are not able to sample directly from Eq. 共A2兲, so that it is necessary to resort to a mixed scheme,17 that extracts samples for ␪ from its full conditional distribution Eq. 共A1兲 and for ␴ 2 from the following approximation of Eq. 共A2兲:

冊

N 共 ␾ ⫺U ␪ 兲 T 共 ␾ ⫺U ␪ 兲 . , ␥ 2⫹ 2 2

In order to ensure that the generated chain converges in distribution to the desired posterior distribution, the samples for ␴ 2 have to be accepted with a suitable probabilistic rule 共for more details see Ref. 17兲. Model 4 For this model we used the Gibbs sampling that requires us to draw samples from the following full conditional distributions:

兺i

˜ ⫺U ˜ ˜␪ 兲共 ␾ ˜ ⫺U ˜ ˜␪ 兲 T 共␾ i i i i

冊

⫺1

.

REFERENCES

⫻exp共 ⫺ ␥ 2 ␴ ⫺2 兲 exp共 ⫺1/2共 ␾ ⫺U T ␪ 兲 T

冉

兺i

⫺1

␳ 1 ⫽ ␳ ⫹N,

⬀Det共 ⌿⫹ ␴ 2 I 兲 ⫺1/2共 ␴ ⫺2 兲 ␥ 1 ⫺1

⌫ ␥ 1⫹

冉 冉

˜ ⫺1˜␪ ⫹ ⫻ ⌺ 0 0

p 共 ␪ 兩 ␴ 2 , ␾ 兲 ⬀p 共 ␪ , ␾ 兩 ␴ 2 兲 ⫽p 共 ␪ 兲 p 共 ␾ 兩 ␪ , ␴ 2 兲 ⫽N 共 ␮ ,⍀ 兲 , 共A1兲 p

i

冊

1

Bailar, J. C. The promise and problems of meta-analysis. N. Engl. J. Med. 337:1–3, 1997. 2 Beal, S. L., and L. B. Sheiner. Nonmem Users Guide. San Francisco, CA: University of California Press, 1998. 3 Byrne, M. M., J. Sturis, R. J. Sobel, and K. S. Polonsky. Elevated plasma glucose 2h postchallenge predicts defects in beta-cell function. Am. J. Physiol. 270:E572–E579, 1996. 4 Carlin, B. P., and T. A. Louis. Bayes and Empirical Bayes Methods for Data Analysis. London: Chapman & Hall, 1996. 5 Cavaghan, M. K., D. A. Ehrmann, M. M. Byrne, and K. S. Polonsky. Treatment with the oral antidiabetic agent troglitazone improves beta-cell responses to glucose in subjects with impaired glucose tolerance. J. Clin. Invest. 100:530–537, 1997. 6 Davidian, M., and D. M. Giltinan. Nonlinear Models for Repeated Measurement Data. New York: Chapman & Hall, 1995. 7 Eaton, R. P., R. C. Allen, D. S. Schade, K. M. Erickson, and J. Standefer. Prehepatic insulin production in man: Kinetic analysis using peripheral connecting peptide behaviour. J. Clin. Endocrinol. Metab. 51:520–528, 1980. 8 Faber, O. K., C. Hagen, C. Binder, J. Markussen, V. K. Nahitani, P. M. Blix, H. Kuzuya, D. L. Horwitz, A. H. Rubenstein, and N. Rossing. Kinetics of human connecting peptide in normal and diabetic subjects. J. Clin. Invest. 62:197–202, 1979. 9 Ferrannini, E., and C. Cobelli. The kinetics of insulin in man. I. General aspects. Diabetes Metab. Rev. 3:335–363, 1987. 10 Gelman, A., F. Bois, and J. Jiang. Physiological pharmacokinetic analysis using population modelling and informative prior distributions. J. Am. Stat. Assoc. 91:1400–1412, 1996. 11 Gelman, A., J. B. Carlin, H. S. Stern, and D. B. Rubin. Bayesian Data Analysis. London: Chapman & Hall, 1995. 12 Geman, S., and D. Geman. Stochastic relaxation, Gibbs distributions, and the bayesian restoration of images. IEEE Trans. Pattern Anal. Mach. Intell. 6:721–741, 1984. 13 Gilks, W. R., S. Richardson, and D. J. Spiegelhalter. Markov

Bayesian Identification of a Population Model of CP Kinetics Chain Monte Carlo in Practice. London: Chapman & Hall, 1996. 14 Hastings, W. K.. Monte Carlo sampling methods using Markov Chain and their applications. Biometrika 57:97–109, 1970. 15 Hovorka, R., E. Koukkou, D. Southerden, J. K. Powrie, and M. A. Young. Measuring insulin secretion using a population model of C-peptide kinetics: Accuracy and required sampling schedule. Diabetologia 41:548–554, 1998. 16 Hovorka, R., M. A. Young, and P. A. Soons. ISEC: A program to calculate insulin secretion. Comput. Methods Programs Biomed. 50:253–264, 1996. 17 Magni, P., R. Bellazzi, and G. De Nicolao. Bayesian function learning using MCMC methods. IEEE Trans. Pattern Anal. Mach. Intell. 20:1319–1331, 1998. 18 Mentre, F., and A. Mallet. Handling covariates in population pharmacokinetics. Int. J. Bio-med. Comput. 36:25–33, 1994. 19 Metropolis, N., A. W. Rosenbluth, M. N. Rosenbluth, M. N. Teller, and A. H. Teller. Equations of state calculations by fast computing machine. J. Chem. Phys. 21:1087–1091, 1953. 20 O’Meara, N. M., J. Sturis, E. Van Cauter, and K. S. Polonsky. Lack of control by glucose of ultradian insulin secretory oscillations in impaired glucose tolerance and in non-insulindependent diabetes mellitus. J. Clin. Invest. 92:262–271, 1993.

21

823

Polonsky, K. S., J. Licinio-Paixao, B. D. Given, W. Pugh, P. Rue, J. Galloway, T. Karrison, and B. Frank. Use of biosynthetic human C-peptide in the measurement of insulin secretion rates in normal volunteers and type i diabetic patients. J. Clin. Invest. 77:98–105, 1986. 22 Shapiro, E. T., H. Tillil, A. H. Rubenstein, and K. S. Polonsky. Peripheral insulin parallels changes in insulin secretion more closely than C-peptide after bolus intravenous glucose administration.J. Clin. End. Metab. 67:1094–1099, 1988. 23 Sparacino, G., and C. Cobelli. A stochastic deconvolution method to reconstruct insulin secretion rate after a glucose stimulus. IEEE Trans. Biomed. Eng. 43:512–529, 1996. 24 Toffolo, G., F. De Grandi, and C. Cobelli. Estimation of beta cell sensitivity from IVGTT C-peptide data. Knowledge of the kinetics avoids errors in modeling the secretion. Diabetes 44:845–854, 1995. 25 Van Cauter, E., F. Mestrez, J. Sturis, and K. S. Polonsky. Estimation of insulin secretion rates from C-peptide levels. Diabetes 41:368–377, 1992. 26 Vonesh, E. F.. Nonlinear models for the analysis of longitudinal data. Stat. Med. 11:1929–1954, 1992. 27 Wakefield, J., and J. Bennett. The bayesian modelling of covariates for population pharmcokinetic models. J. Am. Stat. Assoc. 91:917–927, 1996.