Estimating postprandial glucose fluxes using hierarchical Bayes modelling

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 102–112

journal homepage: www.intl.elsevierhealth.com/journals/cmpb

Estimating postprandial glucose fluxes using hierarchical Bayes modelling Ahmad Haidar a,d , Elizabeth Potocka b , Benoit Boulet d , A. Margot Umpleby c , Roman Hovorka a,∗ a

University of Cambridge Metabolic Research Laboratories, Cambridge, UK MannKind Corporation, Paramus, NJ, USA c Diabetes and Metabolic Medicine, University of Surrey, Guilford, UK d Centre for Intelligent Machines, McGill University, Canada b

a r t i c l e

i n f o

a b s t r a c t

Article history:

A new stochastic computational method was developed to estimate the endogenous glucose

Received 16 June 2011

production, the meal-related glucose appearance rate (Ra meal ), and the glucose disposal (Rd )

Received in revised form

during the meal tolerance test. A prior probability distribution was adopted which assumes

22 December 2011

smooth glucose fluxes with individualized smoothness level within the context of a Bayes

Accepted 25 January 2012

hierarchical model. The new method was contrasted with the maximum likelihood method using data collected in 18 subjects with type 2 diabetes who ingested a mixed meal con-

Keywords:

taining [U-13 C]glucose. Primed [6,6-2 H2 ]glucose was infused in a manner that mimicked

Meal tolerance test

the expected endogenous glucose production. The mean endogenous glucose production,

Bayesian estimation

Ra meal , and Rd calculated by the new method and maximum likelihood method were nearly

Glucose flux

identical. However, the maximum likelihood gave constant, nonphysiological postprandial

Glucose tracer

endogenous glucose production in two subjects whilst the new method gave plausible esti-

Mathematical modelling

mates of endogenous glucose production in all subjects. Additionally, the two methods were compared using a simulated triple-tracer experiment in 12 virtual subjects. The accuracy of the estimates of the endogenous glucose production and Ra meal profiles was similar [root mean square error (RMSE) 1.0 ± 0.3 vs. 1.4 ± 0.7 ␮mol/kg/min for EGP and 2.6 ± 1.0 vs. 2.9 ± 0.9 ␮mol/kg/min for Ra meal ; new method vs. maximum likelihood method; P = NS, paired t-test]. The accuracy of Rd estimates was significantly increased by the new method (RMSE 5.3 ± 1.9 vs. 4.2 ± 1.3; new method vs. ML method; P < 0.01, paired t-test). We conclude that the new method increases plausibility of the endogenous glucose production and improves accuracy of glucose disposal compared to the maximum likelihood method. © 2012 Elsevier Ireland Ltd. All rights reserved.

1.

Introduction

Measurement of postprandial glucose turnover facilitates understanding of physiology in health and under pathophysiological conditions [1–4]. However, accurate measurements are hindered by the nonsteady state of the glucoregulatory

system in the postprandial conditions and ill-conditioning of the inverse problem [5]. The dual tracer dilution methodology, known as the double tracer (DT) technique, has been proposed to asses postprandial glucose fluxes [6] and used by different investigators [3,4,7,8]. Utilizing this technique, one glucose tracer is infused

∗ Corresponding author at: University of Cambridge Metabolic Research Laboratories, Institute of Metabolic Science, Box 289, Addenbrooke’s Hospital, Hills Road, Cambridge CB2 0QQ, UK. Tel.: +44 1223 762 862; fax: +44 1223 330 598. E-mail address: [email protected] (R. Hovorka). 0169-2607/$ – see front matter © 2012 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.cmpb.2012.01.010

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 102–112

intravenously to assess endogenous glucose production (EGP) whilst a second glucose tracer is mixed with the meal to assess the rate of appearance of oral glucose [9]. Infusing the intravenous tracer in a manner that mimics the expected EGP improves the estimation of EGP by reducing the model misspecification error [10,11]. This infusion pattern minimizes the change over time in the specific activity or the tracer-to-tracee ratio (TTR). The triple tracer (TT) technique utilizes a third tracer which is infused in a manner that mimics the expected glucose appearance from the meal (Ra meal ) aiming to improve the estimate of Ra meal [12]. The DT and TT experimental approaches need to be accompanied by an appropriate computational method. The one compartment Steele’s model [6] and the two compartment Radziuk/Mari model [4,11] are most widely used to represent the glucose kinetics when estimating the postprandial glucose turnover. The Radziuk/Mari model has been shown to outperform the Steele model but does not entirely avoid the model misspecification error [12]. To address ill-conditioning, two computational methods have been proposed to be used in connection with the model of glucose kinetics. The traditional method utilizes transformations of model equations [6,11,13] whilst a recent method based on the maximum likelihood (ML) theory uses untransformed model equations [5]. The ML method has been shown to outperform the traditional method [5]. The ML method uses a regularizing technique to smooth glucose fluxes rather than using smoothing glucose measurements by the traditional method. The ML method avoids transformation of the measurement error that may skew estimates of the glucose fluxes with the traditional method [5]. However, when the extent of smoothing is not individualized, the ML method may provide non-physiological constant postprandial EGP in about 10% of cases (observed in our unpublished data). This may limit its utility. In the present work, a new computational method based on the hierarchical Bayes modelling is described. The method assumes a population-based prior to estimate the individualized smoothness level. A clinical study utilizing the DT design is used to contrast the new method with the ML method. Additionally, a comparison of the new method and the ML method is made using simulated data with a TT study design.

2.

Materials and methods

2.1.

Introduction

We aim to estimate (1) the endogenous glucose production, (2) the rate of glucose appearance from the meal, and (3) the rate of disposal of glucose, given unlabelled and labelled glucose concentration measurements. Where possible, time dependency and vector dimensions are dropped for the sake of notational simplicity.

2.2.

Model setup with triple-tracer study design

Without loss of generality, the Radziuk/Mari model for a TT study design is used to illustrate the principles. The two compartment Radziuk/Mari model assumes an irreversible loss

103

from the accessible compartment, a glucose appearance rate in the accessible compartment and glucose exchange between the accessible and nonaccessible compartments. The kinetics of glucose species S, S = E, EM, M and MM are described by the following differential equations: dQ1,S (t) = us (t) − (k01,s (t) + k21 )Q1,S (t) + k12 (t)Q2,S (t) dt dQ2,S (t) = k21 Q1,S (t) − k12 (t)Q2,S (t) dt Q1,S (t) GS (t) = V

(1)

where Q1,S (t) and Q2,S (t) are glucose masses in the accessible and non-accessible compartments, respectively, GS (t) is the plasma glucose concentration, V is the distribution volume of glucose in the accessible compartment, uS (t) is the glucose appearance rate, k12 and k21 are fractional turnover rates. The fractional clearance k01,S (t) is defined as:

k01,S (t) =

⎧ ⎪ ⎨ RdM (t)

if S = M, MM

⎪ ⎩

if S = E, EM

Q1,M (t) RdE (t) Q1,E (t)

(2)

where RdM (t) is the rate of disposal of meal glucose (i.e., the proportion of plasma glucose originating from the meal) and RdE (t) is the rate of disposal of endogenous glucose (i.e., the proportion of plasma glucose originating from endogenous sources). The total unlabelled glucose disposal Rd (t) is equal to summation of RdM (t) and RdE (t). The glucose appearance rate is defined as follows:

⎧ R ⎪ ⎨ a

meal (t)

EGP(t) uS (t) = ⎪ ⎩ uEM (t) uMM (t)

if S = M if S = E if S = EM if S = MM

(3)

where uEM (t) and uMM (t) are the known rates of the intravenously infused tracers.

2.3.

Maximum likelihood method

The principles of the maximum likelihood [5] method are illustrated using a simple exposition to limit notational complexity. Let n be the number of samples. The unknown glucose flux F is represented by a piecewise linear function with F(t) = Fj at the sampling points, where j = 1, . . ., n. Fj s are estimated by maximizing the likelihood LLM LLM = LLMfit + LLMreg

(4)

LLMfit corresponds to the goodness of fit and evaluates the difˆ i and model fit G ¯ i, ference between glucose measurements G the latter obtained from Fj by, for example, the Radziuk/Mari model



LLMfit ∝ −

i

¯i −G ˆ i) (G

2

(5)

104

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 102–112

LLMreg evaluates regularity of the glucose flux Fj



LLMreg ∝ −

ln(ω2 ) +

j

(Fj − Fj−1 )

2

ω2

(6)

After suitable algebraic manipulations, the regularization parameter ω2 can be eliminated and Eq. (7) can be rewritten as

⎛



LLMreg ∝ −ln ⎝

⎞ (Fj − Fj−1 )

2

⎠

(7)

j

where the regularizing parameter ω2 is fully determined by the standard deviations of the measurement error of the labelled and unlabelled glucose concentration measurements. Normally, LLM has one improper and one proper (local) maximum, the former arising when all fluxes Fj are identical and the logarithm term in Eq. (7) reaches infinity; and the latter when Fi values are not identical. In certain circumstances, particularly when Fi with relatively small variations fit well the measured glucose, the proper solution does not exist and only the improper solution is obtained. This results in an unsatisfactory “constant” glucose flux and particularly EGP is affected. The conditions governing the existence of the proper solution are unknown to the authors but it was observed that small changes in glucose measurements may lead to the presence/abolition of the proper solution. For detailed explanation of the maximum likelihood method, we refer the reader to Ref. [5].

2.4.

New method based on Bayesian analysis

The new proposed method utilizes the Bayes theory [14] to calculate the posterior density of unknown fluxes f(u|y,v) defined as: f (u|y, v) =

f (y, v|u)f (u) f (y, v)

functions for a DT study design and/or adopting the Steele model can be obtained using a similar approach. Let n be the number of subjects and m be the number of samples taken in each subject. The prior distribution of the unknown glucose flux f, f = EGP, Ra meal , RdE and RdM , is defined as random walk: fj (i + 1) = fj (i) + RWj (i) RWj (i)∼N(0, εij )

where j = 1, . . ., n, i = 0, . . ., m − 1, N(a,b) is the normal distribution with mean a and precision (precision is defined as an inverse of variance) b. Precisions of the random walks determine the extent of smoothness. The random walks share, in a stochastic manner, the precision. The individual precisions (inverse of variances) are sampled from a log-normal distribution. The definition of εij , corrected for non-equidistant sampling schedules, is as follows: εij =

εj

ti+1 − ti exp(εj )∼N(¯ε, )

(10)

ε¯ ∼ (0.001, 0.001) ∼ (0.001, 0.001) where ti , i = 0, . . ., m − 1 are the sampling times and  (a,b) is the gamma distribution with parameters a and b, ε¯ is the mean of the population distribution and ¯ is the precision of the population distribution. The values of the fluxes at time zero are defined, assuming steady state conditions in glucose appearance and disposal and realizing that meal-related glucose appearance commence with meal ingestion at time zero: j

Ra (0) = 0 EGPj (0)∼N(EGP, 2 ) 2 ∼ (0.001, 0.001)

(8)

(9)

EGP∼N(10, 1−10 )

(11)

j

where f(u) is the prior density function, f(y|v,u) is the likelihood function and f(y,v) is the joint probability acting as a normalizing factor. The posterior distributions contain information about the mean and median as well as associated uncertainty. Measures of uncertainty such as the 75% and 95% credible intervals can be extracted from the posterior distributions. The choice of the prior density f(u) allows incorporation of prior knowledge regarding unknown parameters, e.g. knowledge from the literature or preliminary data. Importantly, the prior can impose regularization to address the ill-conditioning of the inverse problem [15]. The regularizing prior assumes smooth glucose fluxes with individualized smoothness levels and adopts a population-based approach [16].

2.5.

Prior and likelihood functions

In this section, the prior and the likelihood functions that are used to estimate fluxes during TT study adopting the Radziuk/Mari model are described. The prior and the likelihood

RdE (0) = EGPj (0) j

RdM (0) = 0 where EGP and  2 are the mean and the precision of the population distribution of the premeal EGP. The initial conditions of glucose masses of different species are defined assuming steady state conditions: QM1 (0) = QM2 (0) = QMM1 (0) = QMM2 (0) = 0 uEM (0)QE1 (0) RdE (0) QEM1 (0)k21 QEM2 (0) = k12 QE1 (0)k21 QE2 (0) = k12 QE1 (0)∼N(Q¯ E1 , 3 ) QEM1 (0) =

Q¯ E1 ∼N(600, 1−10 ) 3 ∼ (0.001, 0.001)

(12)

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 102–112

105

where QE1 (0), QE2 (0), QM1 (0), QM2 (0), QMM1 (0), QMM2 (0), QEM1 (0) and QEM1 (0) are glucose masses of different glucose species as defined previously, QE1 and  3 are the mean and the precision of the population distribution of the initial endogenous glucose. The time course of glucose masses of different species, i.e. ˆ S (i) as a solution of Eq. (1), is related to the measurements G follows:

4. Solve numerically the differential equations of the Radziuk/Mari model with accompanying initial conditions and individual glucose fluxes as defined in the proposed state ¯ (i+1) . 5. The proposed state ¯ (i+1) is accepted, i.e.  (i+1) = ¯ (i+1) , if ˛, drawn from the standard uniform distribution, satisfies

ˆ EM (i) = QE1 (i) + eEM (i) i = 0, . . . , m G V ˆ M (i) = QM1 (i) + eM (i) i = 1, . . . , m G V ˆ MM (i) = QMM1 (i) + eMM (i) i = 1, . . . , m G V ˆ T (i) = QT1 (i) + et (i) i = 0, . . . , m G V  QT1 (i) = QEM1 (i) + Qe1 (i) + QMM1 (i) + QM1 (i) 1 +

where f is known as a proposal or jumping density [17–19]. Note that P(Y| (i+1) ) and P(Y|¯ (i+1) ) are calculated from the Bayesian likelihood function and P(¯ (i+1) ) and P( (i + 1) ) are calculated from the Bayesian prior distribution, as defined in Section 2.5. If the proposed new state is rejected, then the next state is the same as the current state, i.e.  (i + 1) =  (i) . 6. Go to Step 3.

(13)

1 TTRmeal



where eEM (i), eM (i), eMM (i) and eT (i) are the measurement errors and TTRmeal is the tracer-to-tracee ratio of the tracer in the meal.

2.6.

Implementation details

The new method utilized a Markov chain Monte Carlo (MCMC) approach [17] to obtain the posterior distributions. The MCMC was implemented using WinBUGS version 1.4 [18,19]. The intrinsic Gaussian Conditional Autoregressive (CAR) model was used to implement the random walks (GeoBUGS, Version 1.2, Imperial College, London, UK). The CAR model implements a sum-to-zero constraint on the random effects, which significantly improves mixing of the Markov chain. Differential equations in Eq. (1) were implemented using the WBDiff interface version 1.9.4 (www.winbugs-development.org.uk). Excerpt of the WinBUGS code is given in Appendix B. MCMC methods are a class of algorithms that allow sampling from complex intractable posterior Bayesian distributions. The MCMC methods create a Markov chain that, under fairly general conditions, converges to the posterior distribution [17]. The Metropolis–Hasting (MH) algorithm is a generally applicable scheme to generate the Markov chain although faster approaches such as Gibbs sampler exist [17]. Basic steps of the MH algorithm are given below. Y = [Y(0), . . . , Y(m − 1)], Y(i) = Let ˆ T (i) G ˆ EM (i) G ˆ M (i) G ˆ MM (i)], represents the measured unla[G belled and labelled glucose concentrations, U represents the individual parameters, including RW(i), ε, EGP(0) and QE1 (0), and P represents the set of population parameters, i.e. EGP, Q¯ E1 , ε¯ , ,  2 and  3 . The MCMC is used to obtain the joint posterior distribution of all unobserved stochastic variables U and P conditional on the observed data Y, i.e. p(U,P|Y). 1. Define a Markov chain  (i) , i = 1, 2, . . ., as  (i) = [U(i) P(i) ]. 2. Initialize the Markov chain  (1) by setting suitable initial values for each of the stochastic variables, e.g. glucose fluxes, smoothing parameters, pre-meal individual EGP values, population means, etc. 3. Propose a new state of the Markov chain ¯ (i+1) that depends only on the previous state  (i) .

˛<

P(Y|¯ (i+1) )P(¯ (i+1) )f ( (i+1) |¯ (i+1) ) P(¯ (i+1) |Y)f ( (i+1) |¯ (i+1) ) = P( (i+1) |Y)f (¯ (i+1) | (i+1) ) P(Y| (i+1) )P( (i+1) )f (¯ (i+1) | (i+1) )

After a sufficient number of samples, the Markov chain converges to the joint posterior distribution of the stochastic parameters and all subsequent samples are considered as posterior realizations of the distribution. Empirical statistics of these samples are then used to infer the true values of the unknown parameters.

2.7.

Experimental study

Fig. 1 shows the design of the clinical study. Eighteen subjects with type 2 diabetes treated by insulin (15 males, 3 females, age 55 ± 7 years, BMI 30 ± 3 kg/m2 ) participated. The study was approved by the local Ethics Committee and subjects provided written informed consent. Subjects received a primed, constant, continuous intravenous infusion of [6,6-2 H2 ]glucose (6 mg/kg prime; 0.06 mg/kg/min continuous infusion; Cambridge Isotope Laboratories, Andover, MA, USA) for 7 h prior to the digestion of a liquid mixed meal (540 kcal, 67.5 g carbohydrate, 21 g protein, and 21 g fat) containing 1 g of [U-13 C]glucose (Cambridge Isotope Laboratories). Variable, intravenous insulin (Actrapid, NovoNordisk, Mainz, Germany) and variable 20%dextrose enriched with [6,62 H ]glucose (8 mg/g) were infused over 5 h prior to the meal 2

Fig. 1 – The experimental study protocol utilizing the two-tracer technique and involving ingestion of a labelled liquid mixed meal in subjects with type 2 diabetes.

106

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 102–112

premeal rate commencing at 0, 10, 20, 30, 220, 240, 260, 280, 300, and 400 min relative to the start of the meal. Arterialized venous blood samples were taken at −20, −15, −10, −5, 0, 10, 20, 30, 40, 50, 60, 75, 90, 105, 120, 135, 150, 180, 210, 240, 270, 300, 330, 360, 420 and 480 min, analysed for glucose (Super GL; Hitado Dellecke-Möhnesee, Germany), and derivatized to form the methoxime-trimethylsilyl derivative for gas chromatography mass spectrometry analysis (Agilent 5975C mass-selective detector; Agilent, Woking, UK) to measure ions m/z 319.2 (M + 0), 321.2 (M + 2), 323.2 (M + 4) for the determination of TTR of [6,6-2 H2 ]glucose and [U-13 C]glucose as described previously [20,21]. EGP, Ra meal and Rd profiles were calculated using the ML method as described in Ref. [5]. EGP, Ra meal and Rd profiles were also calculated using the new method employing the Radziuk/Mari model. The measurement errors used in the new method associated with the total glucose and [6,6-2 H2 ]glucose were assumed to be multiplicative with a coefficient of variation (CV) of 3% and 5%, respectively. The measurement error associated with [U-13 C]glucose was assumed to be multiplicative with a 7% CV if tracer concentration was greater than 0.01 mmol/l and additive with a standard deviation of 0.007 mmol/l otherwise. The new method was implemented as described in Section 2.6.

2.8.

Fig. 2 – The concentration of the plasma glucose (top), EGP-mimicking tracer Gem ([6,6-2 H2 ]glucose; middle), and the tracer included in the meal ([U13 C]glucose; bottom) in an experimental study involving the administration of two stable label glucose tracers during ingestion of a liquid mixed meal in eighteen subjects with type 2 diabetes [mean (SD)].

ingestion to maintain a target glucose level of 5.0 mmol/l. The insulin infusion was fixed 90 min prior to the meal ingestion at a rate to maintain normoglycaemia and continued until the end of the study. Insulin lispro bolus (12 U; Eli Lilly, Indianapolis, IN, USA) was administered subcutaneously immediately prior to meal ingestion. Following meal digestion, the rate of infusion of [6,6-2 H2 ]glucose was altered in a manner that mimicked the expected endogenous glucose production; specifically [6,62 H ]glucose was infused in a piecewise constant manner at 2 70%, 60%, 50%, 35%, 45%, 50%, 55%, 65%, and 70% of the

Simulation study

A simulation study was performed to compare the ML method and the new method employing synthetic data generated with the use of a glucoregulatory model, which comprises a submodel of gut absorption, a submodel of insulin secretion and kinetics, and a submodel of glucose kinetics and insulin action. The submodel of the gut absorption adopts a two compartment structure as described by Hovorka et al. [22]. The submodel of insulin secretion assumes a linear relationship between plasma glucose and insulin secretion [2]. A one compartment model of the insulin kinetics is adopted [22]. A submodel of the glucose kinetics assumes a two compartment structure with three insulin action compartments representing insulin action on the glucose distribution/transport, disposal, and EGP suppression [21]. Twelve synthetic subjects underwent a simulated TT experiment with a 75-g carbohydrate meal. Model parameters of the synthetic subjects were sampled from informed distributions [2,22]. The endogenous mimicking (EM) and meal mimicking (MM) tracers were given as piecewise constant infusions using profiles similar to those described by Basu et al. [12]. The simulated experiments resulted in peak gut absorption of 29 ± 6 ␮mol/kg/min (mean ± SD) at 77 ± 20 min postmeal, a maximum suppression of the EGP of 84 ± 12% at 110 ± 57 min and peak glucose disposal of 37 ± 8 ␮mol/kg/min at 88 ± 44 min. The model-derived total plasma glucose concentration and the TTR of the EM, MM, and meal (M) tracers were sampled at −20, −10, 0, 5, 10, 20, 30, 40, 50, 60, 70, 90, 120, 150, 180, 210, 240, 270, 300, 330 and 360 min. A multiplicative uncorrelated random measurement error of 1.5% CV was added to the measurements of total glucose concentration. An additive uncorrelated random measurement error

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 102–112

107

of 1.32 × 10−5 , 2.40 × 10−6 , and 9.68 × 10−4 SD (unitless) was added to the measurements of the TTR of EM, MM, and M tracers, respectively. The random measurement errors were sampled from normal distributions with means defined by the model-derived glucose concentrations. The synthetic profiles of the unlabelled and labelled glucose concentrations were analysed with the maximum likelihood and new methods to reconstruct EGP, Ra meal and Rd profiles. Differences between the actual and reconstructed EGP, Ra meal and Rd profiles were summarized using the root mean square error (RMSE). A paired t-test was used to contrast RMSE calculated with the ML method and the new method. P < 0.05 was considered statistically significant.

3.

Results

3.1.

Experimental study

The concentrations of the plasma glucose and the two glucose tracers are shown in Fig. 2. The basal plasma glucose concentration averaged 4.9 ± 0.27 mmol/l and peaked at 7.1 ± 1.46 mmol/l at 105 min postmeal. The plasma glucose concentration returned to basal values by ∼300 min. The concentration of [6,6-2 H2 ]glucose averaged 0.18 ± 0.03 mmol/l before meal ingestion and declined to a minimum value of 0.067 ± 0.02 mmol/l at 240 min postmeal, and then increased slightly by the end of the experiment. The concentration of [U-13 C]glucose reached peak of 0.066 ± 0.02 mmol/l at 160 min postmeal. The concentration ratio GEM over GE is shown in Fig. 3. Due to the infusion pattern of EM-tracer mimicking expected EGP, the change in the plasma [6,6-2 H2 ]glucose to endogenous glucose concentration ratios were minimized, therefore minimizing non-steady state model-dependent error in estimating EGP [10,11]. The mean EGP, Ra meal and Rd from the Bayesian analysis are shown in Fig. 4. The corresponding mean fluxes obtained with the use of the ML method are also shown in Fig. 4. The experimental study indicates a nearly identical mean EGP, Ra meal , and Rd obtained with the two computational approaches. Unlike the new method, the ML method gave nonphysiological EGP estimates (constant EGP) in two out of the eighteen subjects when improper maxima arose as solutions to the

Fig. 4 – The mean EGP (top panel), Ra meal (middle panel), and Rd (bottom panel) estimated with the new method based on Bayesian technique compared to the estimates obtained using maximum likelihood in an experimental study [mean (SD)].

Fig. 3 – The tracer-to-tracee ratio of Gem /Gt ([6,6-2 H2 ]glucose over [U13 C]glucose) in an experimental study.

maximum likelihood optimization problems. The EGP estimate in Subject 2 calculated using the new method and the estimate obtained using the ML method is shown in Fig. 5. It is worth mentioning that, due to the ill-conditioning, both fluxes

108

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 102–112

Fig. 5 – A sample individual EGP calculated using the method based on Bayesian technique compared to the estimate obtained using maximum likelihood in an experimental study.

fit the data well and within the measurement error. However, with the new method, the population prior pushes the initial EGP value towards the population mean (7.2 ␮mol/kg/min). Similarly, the extent of temporal smoothing of the EGP profiles is individualized and is sampled from a population distribution. Whereas the smoothing parameter ε is infinity in the ML method, it attains a finite value with the new method. Thus the population approach avoided constant non-physiological EGP profiles.

3.2.

Simulation study

The RMSEs associated with the new and ML methods are shown in Table 1. The RMSEs associated with EGP and Ra meal were comparable (P = NS) between the two methods. The RMSE associated with Rd was significantly improved with the new method compared to the ML method (P < 0.01). This shows that the new method is more robust to model misspecification [5]. On visual inspection, the increase in the RMSE with the ML method was attributable to a marked overestimation of the peak value of Rd .

4.

Discussion

A new, fully probabilistic method based on hierarchical Bayes modelling is proposed to estimate glucose fluxes during the meal tolerance test. The new method calculates posterior probability of the fluxes. The posterior probability distributions are information-rich compared to fluxes calculated by deterministic techniques. In addition to mean and median, various credible intervals can be determined. With the double and triple tracer techniques, the glucose disposal Rd appears to be most affected by model misspecification. The new method showed significantly increased accuracy of Rd over the ML method. This can be attributed to

increased robustness to model misspecification due to using the full Bayesian inference with the population-based prior. In about 10% of cases, the ML method gives a constant non-physiological postprandial EGP profile. Constant EGP profiles result from the improper minimum of the optimization problem accompanying the maximum likelihood method. The smoothing parameter defined as the inverse of the variance of the random walks reaches infinity. Even constant profiles are driven by the data and the constant EGP will be attained so as to provide the best possible fit. The fact that nonphysiological constant profiles can still fit the data depicts the extent of ill-conditioning of the computational problem. Relatively large variations in the EGP estimates give small variations in the model fit in comparison to the measurement error. The new method increases robustness of the EGP estimate and avoids a constant postprandial EGP. Individual parameters such as pre-meal EGP value and the variance of the random walk associated with EGP are drawn from population distributions. The parameters of the population distributions have a low precision, and thus the variance of the prior population distribution is large implying vague prior distributions. The posterior distributions of individual parameters such as the pre-meal EGP values and smoothing parameters are driven by the data and the ability of subjects to influence each other through the hierarchical model structure. This is generally known as information sharing among subjects. Despite constant EGP estimates in two subjects by the ML method, the mean EGP profiles did not differ between the two methods. This is because EGP estimates did not differ substantially between the two methods in the remaining sixteen subjects, in concordance with the simulation study, see Table 1. Although EGP estimates were constant in two subjects, the point-wise difference between the two methods in these two subjects was not large enough to induce a substantial difference in the mean profile, see Fig. 5. The hierarchical structure exploits an assumption of an underlying common distribution for individual parameters. In the population setting, individual parameters are estimated simultaneously and are some form of weighted average, combining the data fit and the population parameters [14]. Care needs to be taken when specifying population distributions for a more heterogeneous cohort, e.g. combined analysis of healthy and diabetes subjects. Both methods are based on solid theoretical foundations underpinned by statistical concepts to address ill conditioning of the estimation problem. The new method adds hierarchical structure to allow information sharing among individuals through a common population distribution. Its strength is that it prevents non-physiological constant EGP profiles. Additionally, simulations suggest improved accuracy of Rd and reduced susceptibility to model-misspecification. The drawback is increased computational complexity. Dedicated public-domain packages such as WinBUGS make implementation feasible to wider community. The strength of the ML method is its relative computational simplicity but this is offset by occasional failure to provide physiological EGP profiles. The ML method provides accurate population summaries. On balance, we advocate the use of the new method especially when interest lies in individual EGP estimates.

109

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 102–112

Table 1 – The root mean square error is shown representing the difference between the actual and estimated individual EGP, Ra meal , and Rd during simulated experiments in twelve synthetic subjects. RMSE maximum likelihood method (␮mol/kg/min) RMSE Bayesian method (␮mol/kg/min)

EGP

Ra meal

1.4 ± 0.7 1.0 ± 0.3

2.6 ± 1.0 2.9 ± 0.9

Rd * *

5.3 ± 1.9 4.2 ± 1.3

Values are mean ± SD. P < 0.01, Bayesian method vs. ML method. P = NS, otherwise.

∗

Existing estimation methods use deterministic techniques and assume that model parameters are known with absolute accuracy and that model parameters are the same among subjects. However, some model parameters are not known exactly and different investigators used different values, see for example [5,11,12] where different values of the distribution volume in the Radziuk/Mari model were used. If coherent handling of uncertain model parameters is needed, the new method presents itself as a promising technique as it allows the incorporation of uncertainty in model parameters in a straightforward and theoretically coherent manner. The Steele’s and Radziuk/Mari models have been widely used to calculate postprandial glucose fluxes relying on analytical solutions of system equations. The new method employs a numerical solution to solve the differential equations. This allows alternative models to be explored in a straightforward fashion, i.e. higher order models or models with non-linear structure to reduce the model misspecification error. In the present study, subjects digested liquid mixed meal. Given the common underlying computational and experimental procedures, it is reasonable to assume similar results when digesting solid meals as long as the meal tracer is absorbed in an identical fashion as the meal glucose, i.e. when the principle of tracer indistinguishability holds.

5.

Conclusions

In conclusion, a new computational method to calculate postprandial glucose fluxes has been developed. The method uses Bayesian inference with a prior that assumes smooth glucose fluxes with individualized smoothness levels and adopts a population approach to estimate individual parameters. The new method increases the plausibility of individual EGP

estimates and significantly improves Rd estimates compared to the maximum likelihood method. Ra estimates did not differ between the two methods. We recommend the new method to be used due to its increased robustness to modelmisspecification and the plausibility of its estimates.

Acknowledgements EP was an employee of MannKind Corporation at the time of manuscript preparation and writing. Clinical experiments were performed at Profil, Neuss, Germany. Supported by Juvenile Diabetes Research Foundation, Diabetes UK, and NIHR Cambridge Biomedical Research Centre.

Appendix A. Nomenclature

EGP GS Q1,S

endogenous glucose production (␮mol/kg/min) plasma concentration of glucose species S (mmol/l) amount of glucose species S in the accessible compartment (␮mol/kg) amount of glucose species S in the non-accessible Q2,S compartment (␮mol/kg) glucose appearance rate due to meal ingestion of Ra meal unlabelled glucose (␮mol/kg/min) Rd disposal of unlabelled glucose (␮mol/kg/min) Species E native (unlabelled) glucose originating from the EGP Species EM EGP-mimicking tracer Species M tracer incorporated in the meal Species MM meal-mimicking tracer Species T glucose of all origin (native and labelled) tracer-to-tracee ratio of glucose species S over native TTRS glucose in the sample

110

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 102–112

Appendix B. WinBUGS code The excerpt of WinBUGS code implementing the triple tracer experiment is given below. model { ## definition of priors for (i in 1:M) { ## loop over subjects ## glucose fluxes defined as random walks for (k in 2:N) { ## loop over time ## Ra[i, 1] ... Ra in subject i at the start of the experiment ## Ra.detail[i, k] ... random walk with zero mean using CAR model Ra[i,k] <- Ra.detail[i,k] - Ra.detail[i,1] + Ra[i, 1] Rdm[i,k] <- Rdm.detail[i,k] - Rdm.detail[i,1] + Rdm[i,1] Rde[i,k] <- Rde.detail[i,k] - Rde.detail[i,1] + Rde[i, 1] EGP[i,k] <- EGP.detail[i,k] – EGP.detail[i,1] + EGP[i, 1]

} ## definition of CAR EGP.detail[i,1:N ] ~ Ra.detail[i,1:N ] ~ Rdm.detail[i,1:N ] ~ Rde.detail[i,1:N ] ~

models car.normal(adj[], car.normal(adj[], car.normal(adj[], car.normal(adj[],

weights[], weights[], weights[], weights[],

num[], num[], num[], num[],

RW.EGP.tau[i]) RW.Ra.tau[i]) RW.Rdm.tau[i]) RW.Rde.tau[i])

## initial values for glucose fluxes ## initial EGPs are sampled from a population distribution EGP[i, 1] ~ dnorm(EGP.init.mean,EGP.init.tau) Rde[i, 1] <- EGP[i, 1] Rdm[i, 1] <- 0 Ra[i, 1] <- 0 } ## EGP population mean and precision EGP.init.mean ~ dnorm (10,1.0E-10) EGP.init.tau ~ dgamma(a.tau.gamma, a.tau.gamma) ## population means for the smoothness level log.RW.Ra.tau.tau <- log(RW.Ra.tau.tau) log.RW.Rde.tau.tau <- log(RW.Rde.tau.tau) log.RW.Rdm.tau.tau <- log(RW.Rdm.tau.tau) log.RW.EGP.tau.tau <- log(RW.EGP.tau.tau) RW.EGP.tau.tau ~ dgamma(a.tau.gamma,a.tau.gamma) RW.Ra.tau.tau ~ dgamma(a.tau.gamma,a.tau.gamma) RW.Rdm.tau.tau ~ dgamma(a.tau.gamma,a.tau.gamma) RW.Rde.tau.tau ~ dgamma(a.tau.gamma,a.tau.gamma) ## population precisions for the smoothness levels RW.EGP.tau.tau.tau ~ dgamma(a.tau.gamma,a.tau.gamma) RW.Ra.tau.tau.tau ~ dgamma(a.tau.gamma,a.tau.gamma) RW.Rdm.tau.tau.tau ~ dgamma(a.tau.gamma,a.tau.gamma) RW.Rde.tau.tau.tau ~ dgamma(a.tau.gamma,a.tau.gamma) ## individual smoothness levels sampled from log normal population distributions for (i in 1:M) { RW.Ra.tau[i] ~ dlnorm(log.RW.Ra.tau.tau, RW.Ra.tau.tau.tau) RW.Rde.tau[i] ~ dlnorm(log.RW.Rde.tau.tau, RW.Rde.tau.tau.tau) RW.Rdm.tau[i] ~ dlnorm(log.RW.Rdm.tau.tau, RW.Rdm.tau.tau.tau) RW.EGP.tau[i] ~ dlnorm(log.RW.EGP.tau.tau, RW.EGP.tau.tau.tau) }

## ## ## ## ## ## ##

modeltt (.) is a BlackBox function that solves numerically the differential equations associated with the triple tracer method uem[,] and umm[,] represents the rates of the infused tracers init[,] represents the initial conditions of the differential equations t[] represents the vector of times at which the solution is to be evaluated tol represents the numerical tolerance to be used by the Runge-Kutta algorithm origin represents the time to which the initial condition applies

for(i in 1:M){ ## loop over subjects solution[i, 1:2*N-1, 1:dim]
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 102–112

theta[i, theta[i, theta[i, theta[i, theta[i,

N+ j] 2*N+ j] 3*N+ j] 4*N+ j] 5*N+ j]

<<<<<-

111

Ra[i, j] Rde[i, j] Rdm[i, j] umm[i, j] uem[i, j]

} for (j in 1:2*N-1){ theta[i, 7*N + j]<- t[j] } } ## This section defines the measurement model ## G.meas[,],Gem.meas[,],Gm.meas[,] and Gmm.meas[,] represent the measurements ## Gt.tau[,],Gem.tau[,],Gm.tau[,] and Gmm.tau[,] represent the precisions (1/variance) ## associated with the measurement error ## V is volume of distribution for (i in 1:M){ ## loop over subjects for (j in 1:N){ ## loop over time ## measurements are sampled from normal distribution ## defined by the solution of the differential equations (mean) ## and precision Gem.meas[i,j]~dnorm(Gem[i,j],Gem.tau[i, j]) G.meas[i,j] ~dnorm(Gt[i,j],Gt.tau[i, j]) Gem[i,j] Ge[i,j]

<- solution[i,2*j-1,Qem1L]/V <- solution[i,2*j-1,Qe1L]/V

Gt[i, j]

<- Gm[i,j]+Ge[i,j]+Gmm[i,j]+Gem[i,j]+Gm[i,j]/TTRmeal

} Gm[i, 1] <- 0 Gmm[i, 1]<- 0 for (j in 2:N){ Gm.meas[i,j]~dnorm(Gm[i,j], Gm.tau[i, j]) Gmm.meas[i,j]~dnorm(Gmm[i,j], Gmm.tau[i, j]) Gm[i,j]
solution[i,2*j-1,Qm1L]/V solution[i,2*j-1,Qmm1L]/V

} } }

references

[1] E. Ferrannini, O. Bjorkman, G.A. Reichard Jr., A. Pilo, M. Olsson, J. Wahren, R.A. DeFronzo, The disposal of an oral glucose load in healthy subjects. A quantitative study, Diabetes 34 (1985) 580–588. [2] R. Hovorka, L. Chassin, S.D. Luzio, R. Playle, D.R. Owens, Pancreatic beta-cell responsiveness during meal tolerance test: model assessment in normal subjects and subjects with newly diagnosed noninsulin-dependent diabetes mellitus, The Journal of Clinical Endocrinology and Metabolism 83 (1998) 744–750. [3] A. Thorburn, A. Litchfield, S. Fabris, J. Proietto, Abnormal transient rise in hepatic glucose production after oral glucose in non-insulin-dependent diabetic subjects, Diabetes Research and Clinical Practice 28 (1995) 127–135. [4] J. Radziuk, T.J. McDonald, D. Rubenstein, J. Dupre, Initial splanchnic extraction of ingested glucose in normal man, Metabolism 27 (1978) 657–669. [5] R. Hovorka, H. Jayatillake, E. Rogatsky, V. Tomuta, T. Hovorka, D.T. Stein, Calculating glucose fluxes during meal tolerance test: a new computational approach, American Journal of Physiology-Endocrinology and Metabolism 293 (2007) E610–E619. [6] R. Steele, C. Bjerknes, I. Rathgeb, N. Altszuler, Glucose uptake and production during the oral glucose tolerance test, Diabetes 17 (1968) 415–421. [7] R. Rabasa-Lhoret, Y. Burelle, F. Ducros, J. Bourque, C. Lavoie, D. Massicotte, F. Peronnet, J.L. Chiasson, Use of an alpha-glucosidase inhibitor to maintain glucose

[8]

[9]

[10]

[11]

[12]

[13]

[14] [15]

homoeostasis during postprandial exercise in intensively treated Type 1 diabetic subjects, Diabetes Medicine 18 (2001) 739–744. M.E. Pennant, L.J. Bluck, M.L. Marcovecchio, B. Salgin, R. Hovorka, D.B. Dunger, Insulin administration and rate of glucose appearance in people with type 1 diabetes, Diabetes Care 31 (2008) 2183–2187. R.R. Wolfe, Radioactive and Stable Isotope Tracers in Biomedicine: Principles and Practice of Kinetic Analysis, Wiley-Liss, New York, 1992. C. Cobelli, A. Mari, E. Ferrannini, Non-steady state: error analysis of Steele’s model and developments for glucose kinetics, American Journal of Physiology-Endocrinology and Metabolism 252 (1987) E679–E689. A. Mari, Estimation of the rate of appearance in the non-steady state with a two-compartment model, American Journal of Physiology-Endocrinology and Metabolism 263 (1992) E400–E415. R. Basu, B. Di Camillo, G. Toffolo, A. Basu, P. Shah, A. Vella, R. Rizza, C. Cobelli, Use of a novel triple-tracer approach to assess postprandial glucose metabolism, American Journal of Physiology-Endocrinology and Metabolism 284 (2003) E55–E69. G. DeNicolao, G. Sparacino, C. Cobelli, Nonparametric input estimation in physiological systems: problems, methods, and case studies, Automatica 33 (1997) 851–870. P. Congdon, Bayesian Statistical Modelling, 2nd ed., John Wiley & Sons, Chichester, England/Hoboken, NJ, 2006. A.N. Tikhonov, Numerical Methods for the Solution of Ill-Posed Problems, Kluwer Academic Publishers, Dordrecht/Boston, 1995.

112

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 102–112

[16] L.B. Sheiner, B. Rosenberg, V.V. Marathe, Estimation of population characteristics of pharmacokinetic parameters from routine clinical data, Journal of Pharmacokinetics and Biopharmaceutics 5 (1977) 445–479. [17] W.R. Gilks, S. Richardson, D.J. Spiegelhalter, Markov Chain Monte Carlo in Practice, Chapman & Hall, Boca Raton, FL, 1998. [18] D.J. Lunn, A. Thomas, N. Best, D. Spiegelhalter, WinBUGS—a Bayesian modelling framework: concepts, structure, and extensibility, Statistics and Computing 10 (2000) 325–337. [19] T.A. Spiegelhalter D, N. Best, D. Lunn, WinBugs User Manual, Medical Research Council Biostatistics Unit, Cambridge, UK, 2003. [20] F. ShojaeeMoradie, N.C. Jackson, R.H. Jones, A.I. Mallet, R. Hovorka, A.M. Umpleby, Quantitative measurement of

3-O-methyl-d-glucose by gas chromatography mass spectrometry as a measure of glucose transport in vivo, Journal of Mass Spectrometry 31 (1996) 961–966. [21] R. Hovorka, F. Shojaee-Moradie, P.V. Carroll, L.J. Chassin, I.J. Gowrie, N.C. Jackson, R.S. Tudor, A.M. Umpleby, R.H. Jones, Partitioning glucose distribution/transport, disposal, and endogenous production during IVGTT, American Journal of Physiology-Endocrinology and Metabolism 282 (2002) E992–E1007. [22] R. Hovorka, V. Canonico, L.J. Chassin, U. Haueter, M. Massi-Benedetti, M.O. Federici, T.R. Pieber, H.C. Schaller, L. Schaupp, T. Vering, M.E. Wilinska, Nonlinear model predictive control of glucose concentration in subjects with type 1 diabetes, Physiological Measurement 25 (2004) 905–920.