ASSESSMENTS of postprandial glucose turnover facilitate understanding of carbohydrate physiology in health and diseases (7, 22, 30). In the postprandial state, endogenous glucose production (EGP) and glucose disposal (Rd) vary with time together with the appearance of glucose originating from the meal (Ra meal). Direct measurements of postprandial turnover are generally not possible, and indirect assessments using the tracer dilution methodology are hindered by the nonsteady state of the glucoregulatory system (1) and ill-conditioning (small measurement errors propagate into large estimation errors) of the computational problem (12).

The dual-tracer (DT) technique is the original method to assess postprandial glucose fluxes (29) and has been adopted in a range of studies (19, 21, 22, 30). The DT technique uses one glucose tracer infused intravenously to measure EGP while a second glucose tracer mixed with the meal allows oral glucose appearance to be determined (35). Infusing the intravenous tracer in a manner that mimics the expected EGP pattern minimizes the change over time in the specific activity or the tracer-to-tracee ratio (TTR) and improves accuracy of EGP estimates (3, 18).

The triple-tracer (TT) technique utilizes a third tracer, which is infused in a manner that mimics the expected glucose appearance from the meal (1). The TT technique has been proposed to be the gold standard to measure postprandial glucose appearance (1). The TT technique provides a model-independent estimate of Ra meal if TTR of the meal tracer over the meal-mimicking tracer does not change over time. However, for practical reasons, variations in the TTR are inevitable, and errors of unknown magnitude arise resulting from model misspecification and ill conditioning.

Validation of the TT and DT methods against an independent standard is part of our research to measure gut absorption rates of meals containing complex carbohydrates in subjects with type 1 diabetes (T1D) (6). In the present study, we assessed the accuracy of the TT and DT techniques by contrasting the administered rate of intravenous dextrose against the reconstructed rate. The delivery pattern of intravenous dextrose mimicked a meal-derived glucose appearance serving as an independent standard. We used an advanced computational approach based on population stochastic modeling to reduce the effect of ill conditioning.

From 1800 until 0200, a variable intravenous insulin infusion was given to achieve plasma insulin concentrations comparable to that resulting from a subcutaneous bolus of rapid-acting insulin analog. The total amount of variable insulin corresponded to an insulin bolus associated with a meal containing 120 g of carbohydrates. The calculations utilized individualized insulin-to-carbohydrate ratios. Starting at 1800, a variable intravenous 20% dextrose infusion (the independent standard mimicking appearance of exogenous glucose from a standard meal) enriched with [U-13C]glucose (to represent the tracer that would normally be added to the meal; 13 mg/g, Cambridge Isotope Laboratories) was infused over 8 h to replicate post-meal plasma glucose excursions. Postprandial glucose patterns were taken from a previous study in young subjects with T1D who consumed a standard meal containing 120 g of carbohydrates (6). An adaptive model-predictive controller (gMPC, version 1.0.2; University of Cambridge, Cambridge, UK) was used to adjust the dextrose infusion utilizing information about subject's total daily insulin dose, the predetermined intravenous insulin delivery, and the desirable plasma glucose profile.

From 1000, venous blood samples were taken every 10–15 min for the determination of plasma glucose. Samples were taken at 1500, 1515, and 1530 to determine background glucose enrichments. From 1700, venous blood samples were taken every 10–30 min for the determination of plasma insulin, [U-13C]glucose, [6,6-2H2]glucose, and [U-13C; 1,2,3,4,5,6,6-2H7]glucose. The samples were immediately centrifuged and separated. Plasma glucose was measured immediately while other samples were stored at −80°C until assayed.

TTRs were calculated using a variation of a method described previously (12, 14). In brief, ions m/z M+0, M+2, M+4, and M+8 were used to calculate TTRs of [6,6-2H2]glucose, [U-13C]glucose, and [U-13C; 1,2,3,4,5,6,6-2H7]glucose corrected for recycled glucose using ions m/z M+0, M+3, and M+5. The calculations accounted for spectra overlap (23). It was assumed that [U-13C; 1,2,3,4,5,6,6-2H7]glucose is recycled equally into glucose molecules with m/z of M+4 and M+5 and that [U-13C]glucose is recycled equally into glucose molecules with m/z of M+2 and M+3 (15, 16). A system of algebraic equations was solved analytically, and the solution was reduced and validated using Mathematica (Wolfram Research, Champaign, IL). The endogenous glucose concentration was calculated using a model-independent method as described previously (1).

A stochastic modeling approach was used to reconstruct the exogenous dextrose infusions (Ra exo) and to calculate EGP and Rd. The DT method utilized the plasma concentration of unlabeled glucose, [6,6-2H2]glucose, and [U-13C]glucose. The TT methods utilized additionally the plasma concentration of [U-13C; 1,2,3,4,5,6,6-2H7]glucose. A hierarchical Bayes model was created implementing the Radziuk/Mari two-compartment model of glucose kinetics (18, 22). The method used Bayesian inference with a regularizing prior that assumes smooth glucose fluxes (12), with individualized smoothness levels drawn from a population distribution (27). The method has been described previously (9), and additional details are given in the appendixes.

Statistical analysis.

Values are presented as means ± SD or median (interquartile range) as appropriate. The difference between the actual and reconstructed dextrose infusion rates was determined using the root mean square error (RMSE).

RMSE=∑i=1N(Ra exo,i−Ra exo,i¯)2N

where Ra exo,i¯ is the actual dextrose infusion rate and Ra exo,i is its estimate at time i. The overall dextrose recovery was calculated. The EGP and Rd profiles were summarized using the partial area under the curve (AUC) from time 0 to 60 min and from time 0 to 480 min. A paired t-test was used to contrast differences between the TT and DT methods. P < 0.05 was considered statistically significant.

RESULTS

Plasma glucose and plasma insulin.

Plasma glucose and plasma insulin concentrations are shown in Fig. 2. The plasma glucose concentration before the start of dextrose infusion averaged 5.7 ± 1.1 mmol/l and peaked at 11.3 ± 3.7 mmol/l after 120 min, and then decreased slowly to around 8.0 mmol/l by the end of the experiment. The concentration of plasma insulin was 96 (57–197) pmol/l before the start of the dextrose infusion and peaked at 879 (726–986) pmol/l after 70 min. The plasma insulin concentration returned to basal level by 360 min.

TTRs.

The TTRs are shown in Fig. 3. Due to the infusion pattern of [6,6-2H2]glucose and [U-13C; 1,2,3,4,5,6,6-2H7]glucose, the changes in the ratios [6,6-2H2]glucose over endogenous glucose and [U-13C; 1,2,3,4,5,6,6-2H7]glucose over [U-13C]glucose were reduced, minimizing the model-dependent error when estimating EGP and Ra exo by the TT method.

The ratio [6,6-2H2]glucose over [U-13C]glucose is shown in Fig. 3. The [6,6-2H2]glucose was infused in a manner that mimicked the expected EGP, and this led to marked changes in this particular ratio, potentially leading to a model-dependent error in the estimation of Ra exo by the DT method.

Reconstructing dextrose infusion.

Figure 4 shows a sample dextrose infusion profile and its reconstructed profiles by the TT and DT techniques.

The mean dextrose infusion and its reconstructed patterns are shown in Fig. 5. The calculations indicate nearly identical mean profiles obtained with the two methods, documenting the ability of the two methods to reliably estimate exogenous glucose appearance. The RMSE associated with the mean dextrose infusion profile is shown in Table 1.

Estimation errors over time are shown in Fig. 5. The DT technique tended to underestimate the early portion of the dextrose infusion rate. This could be ameliorated by higher sampling frequency. The estimation error of the DT method tended to be less than zero, in concordance with the observation that the overall dextrose recovery was slightly underestimated.

EGP and glucose disposal.

Figure 6 documents a nearly identical EGP and Rd obtained with the two methods. The plot of differences shown in Fig. 6 suggests that, on an individual basis, EGP estimates by the two methods were nearly identical as documented by the low variability of the differences at each time point. Differences in individual estimates of Rd occurred. The TT technique provided slightly higher estimates of Rd at 30 and 180 min.

Table 2 shows partial AUCs associated with the EGP and Rd profiles. AUCs of EGP were statistically lower with the DT technique, but the difference lacked clinical significance.

Fractional clearance rate.

As expected, the fractional clearance rate obtained by the two methods with the use of [6,6-2H2]glucose were virtually superimposable (Fig. 7). The two methods are fitting the same [6,6-2H2]glucose concentration. Minor differences in the fractional glucose clearance resulted from the way the parameter estimation process was set up utilizing simultaneous fit to all glucose species (see Appendix A). We also confirm that the TT and DT methods were equivalent when estimating EGP. The time course of the factional clearance rate calculated with [U-13C; 1,2,3,4,5,6,6-2H7]glucose by the TT method was slightly higher during the initial study period. This likely reflects a high relative measurement error associated with low concentrations of [U-13C; 1,2,3,4,5,6,6-2H7]glucose at the start of the study. This overestimation is linked to a higher [U-13C]glucose turnover calculated by the TT method contributing to a slightly higher overall Rd (Fig. 6).

DISCUSSION

The present study suggests that the TT and DT techniques provide accurate estimates of exogenous glucose appearance. The mean dextrose infusion was accurately reconstructed. The total amount of delivered dextrose was well recovered by the TT technique, whereas the DT technique underestimated the amount by 8%.

We adopted a unique validation method mimicking postprandial conditions and utilizing an independent standard represented by the intravenous dextrose infusion. A concomitant intravenous insulin infusion achieved typical postprandial insulin excursions in subjects with T1D. We focused on T1D because accurate estimation of gut absorption is required for optimum timing and amount of prandial insulin (2, 25). Optimal prandial insulin dosing is important to minimize postprandial glucose excursions, which are important determinants of overall glucose control in T1D (24). Examining T1D avoids the need to consider endogenous insulin secretion, facilitating a simpler characterization of the relationship between plasma insulin and glucose excursions. In combination with mathematical modeling, estimates of gut absorption can be used to explore, in silico (34), alternative insulin dosing strategies including closed-loop systems (5, 10, 11, 13, 33). Given the common underlying computational principles and experimental procedures, it is reasonable to assume applicability of present results to other populations.

On theoretical grounds, the model misspecification error is the only source for differences between the DT and TT techniques. If the Radziuk/Mari model provided an accurate description of glucose kinetics, the two methods would return identical results. In the early portion of the experiment, when dextrose delivery was rapidly increasing, the DT method underestimated the dextrose infusion. The overall amount of dextrose was slightly underestimated by the DT method, in concordance with results obtained by others (1, 32); otherwise, results provided by the two techniques were comparable. The present study thus suggests that the two-compartment Radziuk/Mari model provides a reasonable approximation of glucose kinetics during meal-like experimental scenarios. This is supported by the observation that the fractional turnover rates were comparable apart from the initial part of the study. We argue that the DT technique is the preferred approach, given the reduced cost and reduced experimental and analytic complexity, unless highly accurate measurements of the early postprandial period are required.

The DT approach is equally suitable to estimate EGP as is the TT approach. From a physiological viewpoint, the results were comparable between the two techniques despite the statistical significance in partial AUCs. Indirectly, the present study also validates estimates of EGP. Given the comparability of the computational approach to estimate EGP and dextrose infusion, it is reasonable to assume validity for estimates of EGP, given that the dextrose infusion was accurately reconstructed. The main difference is that EGP originated endogenously but the dextrose infusion exogenously, but this should not affect the validity of the calculations.

The infusion pattern of EGP-mimicking tracer was designed to minimize changes over time of the TTR GEM/GE but led to considerable fluctuations in the ratio GEM/GM (Fig. 3), possibly accentuating model misspecification errors by the DT method. Despite glucose fluxes calculated by the DT and TT methods being similar, they indicate that the model misspecification error may have been mitigated by the use of the novel computational method. Further studies are warranted to determine whether the amount of exogenous glucose would be less underestimated by the DT method if the infusion of the EGP-mimicking tracer was kept constant, as is the case in many (21, 30), but not all (20), DT studies. We suggest that the infusion pattern of EGP-mimicking tracer should be designed to match the research question. Studies focusing on measuring EGP may adopt a pattern similar to that adopted in the present study.

The low initial concentrations of the meal-mimicking tracer and the meal tracer at the start of the experiment caused inaccurate estimation of the initial portion of the fractional clearance rate k01,MM (Fig. 7). However, unlike the traditional calculation approach using an analytic solution and the fractional clearance rate (1, 18), our approach avoids the use of fractional clearance rate, shown here for illustrative purposes, and utilizes model fit to measured concentrations. This renders the calculations more robust and less prone to ill conditioning. As the measurement errors associated with the meal tracer and the meal-mimicking tracer at the early portion of the experiment are large compared with actual measurements, the initial part of the dextrose infusion estimate will be less dependent on the measured concentrations and will more reflect the assumption of “smooth glucose fluxes”. This is a drawback of the TT study design; the DT approach is not affected by similar considerations.

The novel computational approach may have increased the accuracy of the calculations. The approach assumed that glucose fluxes were smooth to limit the effect of ill conditioning (12, 31). The extent of smoothing was shared among individuals by assuming a population distribution from which individual smoothing estimates were drawn. This “pulled” individual smoothing factors closer to their mean value (27) and avoided aberrations and oscillations, which can be present when data are processed separately. Additionally, this novel approach used an appropriate model of the measurement error. The traditional method (1), although computationally simpler, transforms data, which leads to transformations of the measurement error. This is of particular concern at the initial part of the study when concentrations of meal tracer are low and the associated relative measurement errors are high. This was addressed previously by excluding the initial portion of the data from calculations (1), but all data can be utilized with the present approach.

The computational approach estimates all glucose fluxes simultaneously by fitting measurements of all glucose species in all subjects within a single parameter estimation run (Appendix A). This facilitates interaction between fits to different glucose species. This also explains why the DT and TT methods provide slightly different estimates of endogenous glucose production and fractional clearance rate k01EM while fitting identical EGP-mimicking tracer concentrations. The TT method will additionally be fitting the meal-mimicking tracer concentration, and this will propagate slightly into estimates of EGP and the fractional clearance rate.

In conclusion, our results suggest that triple-tracer and dual-tracer techniques combined with advanced computational methods can measure accurately and reliably glucose appearance during postprandial-like conditions. The triple-tracer technique tends to outperform slightly the dual-tracer technique, but the latter benefits from reduced experimental and analytic complexity.

GRANTS

This work was supported by Juvenile Diabetes Research Foundation (nos. 22-2006-1113, 22-2007-1801, 22-2009-801, 22-2009-802), Diabetes UK (BDA07/0003549), the National Institute of Diabetes and Digestive and Kidney Diseases (1R01 DK-085621), Medical Research Council Centre for Obesity and Related metabolic Diseases, and the National Institute for Health Research Cambridge Biomedical Research Centre.

DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

Angie Watts provided laboratory support. The Diabetes Research Network Laboratory Wales (Dr. Steve Luzio) measured plasma insulin. Tomas Hovorka devised TTR calculations and validated the reduced formulae. We acknowledge support by the staff at the Addenbrooke's Wellcome Trust Clinical Research Facility. We are grateful to study volunteers for their participation.

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.

Prior and Likelihood Functions

The prior and the likelihood functions with the TT method adopting the Radziuk/Mari model are described below. The prior and the likelihood functions for the DT method can be obtained using a similar approach. In what follows, Ra exo represents the dextrose infusion, RdE represents the disposal of endogenous glucose, and RdM represents the disposal of glucose originating from the dextrose infusion.

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 exo, RdE, and RdM, is defined as random walk:

[Formula ID: FDA1](A1)

fj(i+1)=fj(i)+RWj(i)

RWj(i)~N(0,εji)

where j = 1, . . . , n, i = 0, . . . , m − 1, N(a,b) is the normal distribution with mean a and precision (i.e., inverse of the variance) b. The flux f is assumed to be piecewise constant for Ra exo, and piecewise linear for EGP, RdE, and RdM. The random walk enforces temporal smoothing (4) with an extent of smoothing determined by the precision of the random walk εji, corrected for nonequidistant sampling schedules and sampled from a log-normal distribution to assure positivity:

εji=εjti+1−ti

exp(εj)~N(ε,τ)

ε~Γ(0.001,0.001)

τ~Γ(0.001,0.001)

where t1,i=0, . . . , m−1, are the sampling times and Γ(a,b) is the gamma distribution with parameters a and b, ε is the mean precision of the population distribution and τ is the precision of the population distribution representing inter-patient variability. Gamma distributions are typically used to define noninformative priors for precisions (4). A log-normal population distribution is used to ensure positivity of the precisions. The values of the fluxes at time 0 are defined, assuming steady-state conditions of glucose appearance and disposal and noting that exogenous glucose infusion commences at time 0, as:

Ra exoj(0)=0

EGPj(0)~N(EGP,τ2)

τ2~Γ(0.001,0.001)

EGP~N(10,1−10)

RdEj(0)=EGPj(0)

RdMj(0)=0

where EGP and τ2 are the mean and the precision of the normal population distribution of the EGP at time 0. EGP is drawn from a vague normal distribution characterized by a low precision and thus a large variance (variance is an inverse of precision). The initial conditions of glucose masses are defined, assuming steady-state conditions, as follows:

QM1(0)=QM2(0)=QMM1(0)=QMM2(0)=0

QEM1(0)=uEM(0)⋅QE1(0)/RdE(0)

QEM2(0)=QEM1(0)⋅k21/k12

QE2(0)=QE1(0)⋅k21/k12

QE1(0)~N(QE1,τ3)

QE1~N(600,1−10)

τ3~Γ(0.001,0.001)

where QE1(0), QE2(0), QM1(0),QM2(0), QMM1(0), QMM2(0), QEM1(0), and QEM2(0) are glucose masses of different glucose species, QE1 and τ3 are the mean and the precision of the normal population distribution of the initial mass of endogenous glucose and are drawn from vague prior distributions.

The time courses of glucose masses of different species as calculated by the Radziuk/Mari model are related to the measurements of tracer and total plasma glucose concentrations G^EM(i), G^M(i), G^MM(i), and G^T(i) as follows:

G^EM(i)=QEM1(i)/V+eEM(i)

G^M(i)=QM1(i)/V+eM(i)

G^MM(i)=QMM1(i)/V+eMM(i)

G^T(i)=QT1(i)/V+eT(i)

QT1(i)=QEM1(i)+QE1(i)+QMM1(i)+QM1(i)⋅(1+1/TTRdex)

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

The method estimates individual parameters such as vertices of random walks RW(i) for each of the four time-varying profiles EGP, Ra exo, RdE, and RdM, one “smoothness” parameter ε per random walk (all together four smoothness parameters), EGP at time 0 EGP(0) and the mass of glucose in the accessible compartment at time 0 QE1(0), and the population parameters EGP,QE1,ε,, τ, τ2, and τ3. In total there are 147 parameters per individual and 6 population parameters. Individual parameters from all subjects and population parameters are estimated simultaneously during a single parameter estimation run (in total, 2,358 parameters) using all available measurements of the total and labeled glucose concentrations.

Implementation Details

We use the Markov chain Monte Carlo (MCMC) method (8) to obtain the posterior distributions of the fluxes. The MCMC was implemented using WinBUGS version 1.4 (17, 28), extended by the WBDiff package version 1.9.4 (MRC Biostatistics Unit, Cambridge, UK). The latter is used to solve numerically the differential equations defining the Radziuk/Mari model. WinBUGS is a public domain package facilitating Bayesian estimation using MCMC methods. An outline of MCMC principles is provided below.

For the purposes of parameter estimation, measurement errors were assumed to be normally distributed with zero mean. The measurement errors associated with the total glucose and [6,6-2H2]glucose were assumed to be multiplicative with coefficients of variation (CV) of 2 and 5%, respectively. The measurement errors associated with [U-13C]glucose and [U-13C; 1,2,3,4,5,6,6-2H7]glucose were assumed to be multiplicative with 5% CV if tracer concentrations were greater than 0.015 mmol/l and 0.01 mmol/l, respectively, and additive with standard deviations of 0.0075 and 0.005 mmol/l, respectively, otherwise. The additive measurement errors were determined empirically and expressed that, at low tracer concentrations, occurring at the early portion of the experiment after the start of tracer infusions, the measurement instrumentation has precision independent of the measured value and that instrumentation detection limits apply.

APPENDIX B: RADZIUK/MARI MODEL SETUP WITH THE TRIPLE-TRACER METHOD

The kinetics of glucose species S, S = E, EM, M, and MM are described by the following differential equations.

dQ1,S(t)dt=uS(t)−(k01,S(t)+k21)Q1,S(t)+k12Q2,S(t)

dQ2,S(t)dt=k21Q1,S(t)−k12Q2,S(t)

GS(t)=Q1,S(t)V

where Q1,S(t) and Q2,S(t) are glucose masses in the accessible and nonaccessible 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, and k12 and k21 are fractional turnover rates. The parameters of the Radziuk/Mari model were identical across all subjects, and glucose species and were assigned standard values k12 = 0.07 min−1, k21 = 0.05 min−1, and V = 160 ml/kg (12).

The fractional clearance k01,S(t) is defined as:

[Formula ID: FDA2]A2

k01,S(t)={RdM(t)/Q1,MifS=M,MMRdE(t)/Q1,EifS=E,EM

where RdM(t) is the rate of disposal of meal glucose, and RdM(t) is the rate of disposal of endogenous glucose. The use of two fractional clearance rates instead of a single rate is the main feature of the TT method. Whereas the DT method assumes identical k01 across all species, the TT method assumes one k01 to describe kinetics of the M-tracer and the MM-tracer and another k01 to describe kinetics of the EM-tracer and endogenous glucose. The difference between estimates of the two fractional clearance rates by the TT method is an indicator of model misspecification. However, at the early part of the experiment, the fractional clearance rate associated with the M-tracer and the MM-tracer is usually poorly estimated due to low inaccurate tracer concentrations invalidating meaningful comparison of the two fractional rates during this specific period.

The total unlabeled glucose disposal Rd(t) is equal to summation of RdM(t) and RdE(t). The glucose appearance rate is defined as follows:

uS(t)={Ra exo(t)if S=MEGP(t)ifS=EuEM(t)if S=EMuMM(t)if S=MM

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

APPENDIX C: BASIC DESCRIPTION OF MONTE CARLO MARKOV CHAIN METHODS

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 (8). The Metropolis-Hasting algorithm is a generally applicable schema to generate the Markov chain, although faster approaches such as Gibbs sampler exist (8). Basic steps of the Metropolis-Hasting algorithm applied to present computational problem are given below.

Let the vectorY=[Y(0),...,Y(m−1)],Y(i)=[G^T(i)G^EM(i)G^M(i)G^MM(i)]

represent the measurements of total and labeled glucose concentrations, U represents the vector of individual parameters such as RW(i), ε, EGP(0), and QE1(0), and P represents the vector of population parameters EGP,QE1, ε, τ, τ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). The steps are as follows:

where f is known as a proposal or jumping density (8, 17, 28). 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. If the proposed new state is rejected, then the next state is the same as the current state, i.e., θ(i+1) = θ(i).

Go to Step 3.

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. Medians of posterior distributions are then used to infer point estimates of unknown parameters. In the present study, the Fisherian approach was utilized to calculate population mean (or median as appropriate) of glucose fluxes across the individuals.