## Abstract

A new calculation method is proposed to quantify the endogenous glucose production (EGP), the glucose appearance rate due to meal ingestion (R_{a meal}), and the glucose disposal (R_{d}) during a three-tracer study design. The method utilizes the maximum likelihood theory combined with a regularization method to achieve a theoretically coherent computational framework. The method uses the two-compartment formulation of the glucose kinetics. Instead of assuming smoothness of unlabeled and labeled glucose concentrations, the method assumes that the EGP, the R_{a meal}, and the fractional glucose clearance are smooth, increasing plausibility of their individual estimates. The method avoids transformation of the measurement errors, which may skew the estimates of the EGP, R_{a meal}, and R_{d} with the traditional approach. Finally, the sequential nature of the calculations is replaced by calculating the EGP, R_{a meal}, and R_{d} in “one go” to avoid the propagation of the errors from the EGP and R_{a meal} into R_{d}. An example study is shown demonstrating the utility of the approach. A better performance of the new method is demonstrated in a simulation study.

- glucose kinetics
- gut absorption
- endogenous glucose production
- glucose tracer
- mathematical modeling

the postprandial glucose turnover constitutes essential information about physiology in health and pathophysiology of disease state (6, 13, 17). Its assessment is hampered by the nonsteady nature of the glucoregulatory system in the postprandial state. Both plasma glucose and insulin vary considerably.

The assessment calls for the administration of at least two glucose tracers (19). One tracer is infused intravenously to quantify glucose disposition. The second tracer is included in the meal to trace the appearance of oral glucose. Recent work (1, 19) indicates that a three-tracer study design is preferable to minimize the estimation error caused by the misspecification of the model of glucose kinetics. It has been shown (3, 12) that administering glucose tracer in the format mimicking the endogenous glucose production (EGP) to achieve a constant specific activity or tracer-to-tracee ratio (TTR) reduces the model misspecification error when the EGP is calculated. The same principle also applies to the estimation of the glucose appearance from the meal (R_{a meal}), which is mimicked by the third tracer.

Moving on from the one-compartment Steele's model (16) to the Radziuk/Mari two-compartment model (2CM) (12, 13) further improves the accuracy of calculations (1, 19). The main reason for this is the reduction of the model misspecification error. Thus, the three-tracer study design with a 2CM structure appears to be the gold standard method.

Whereas the experimental design has evolved and the model structure has become more complex, the computational method has remained virtually unchanged. It consists of smoothing unlabeled and labeled glucose measurements and inserting them or their ratios into closed-form formulae, which iteratively calculate the EGP, R_{a meal}, and glucose disposal (R_{d}). The calculations are straightforward but are based on assumptions requiring critical appraisal. First, it is assumed that the unlabeled and labeled glucose concentrations are “smooth”. This may be valid for the unlabeled glucose trace, but, because the administration of the labeled glucose in the three-tracer study design is either piecewise constant or piecewise linear, this assumption is difficult to justify for the labeled glucose concentrations. Second, the closed-form formulae transform the measurement error, and this transformation may affect the accuracy of calculation specifically for the meal absorption, given the way the meal-mimicking tracer is delivered. Finally, due to the sequential nature of computation, the R_{d} is calculated from the EGP and R_{a meal}, and thus any error in the EGP and R_{a meal} will propagate into the estimation of the R_{d}. For example, if a temporal underestimation of the R_{a meal} occurs due to a large measurement error at a given time, this error will be compensated by a temporal underestimation of the R_{d} to maintain the mass balance irrespective of whether physiological assumption of R_{d} smoothness is upheld.

In the present work, a new calculation method is described. The method assumes that the underlying metabolic fluxes, the EGP and R_{a meal}, and the glucose clearance are smooth. The method avoids the transformation of the measurement error and avoids the sequential nature of the calculations, making all estimation in “one go” to avoid the propagation of the estimation error from the EGP and the R_{a meal} into the R_{d}. The method can be implemented in a spreadsheet. A sample experimental study demonstrates the utility of the approach and contrasts it with the traditional calculations. Additionally, a comparison of the traditional and the new methods is made, adopting a simulation study.

### Glossary

- D
_{s} - Amount of glucose species S in the meal (g)
- EGP
- Endogenous glucose production (μmol·kg
^{−1}·min^{−1}) - G
_{s} - Plasma concentration of glucose species S (mmol/l)
- Q
_{1,S} - Amount of glucose species S in the accessible compartment (μmol/kg)
- Q
_{2,S} - Amount of glucose species S in the nonaccessible compartment (μmol/kg)
- R
_{a meal} - Glucose appearance rate due to meal ingestion of native glucose (μmol·kg
^{−1}·min^{−1}) - R
_{d} - Disposal of native glucose and all labeled glucose species (μmol·kg
^{−1}·min^{−1}) - Species E
- Native (unlabeled) glucose originating from the EGP
- Species EM
- EGP-mimicking tracer ([1,2,3,4,5,6,6-
^{2}H_{7}]glucose) - Species M
- Tracer incorporated in the meal ([2-
^{2}H_{1}]glucose) - Species MM
- Meal absorption-mimicking tracer ([U-
^{13}C;1,2,3,4,5,6,6-^{2}H_{7}]glucose) - Species N
- Native (unlabeled) glucose incorporated in the meal
- Species T
- Glucose of all origin (native and labeled)
- TTR
_{s} - Tracer-to-tracee ratio of glucose species S over native glucose in the sample
*u*_{s}- Appearance of glucose species S (μmol·kg
^{−1}·min^{−1}) with*u*_{N}= R_{a meal}+ EGP and*u*_{E}= EGP - V
_{1} - Glucose distribution volume in the accessible compartment (ml/kg)

## METHODS

#### Model setup with three-tracer study design.

Following the work by Radziuk et al. (13) and Mari (12), the kinetics of glucose species S, S = E, EM, M, MM, and T, is described by a 2CM, with an irreversible loss from the accessible compartment governed by a time-varying fractional clearance rate, *k*_{01,s}, a time-varying appearance rate, *u*_{s}, and time-invariant fractional turnover rates *k*_{21} and *k*_{12} (wherever possible, the time dependency is dropped for the sake of notational simplicity). (1) (2) (3) where Q_{1,S} and Q_{2,S} are glucose masses in the accessible and nonaccessible compartments, respectively, G_{S} is the plasma glucose concentration, and V_{1} is the time-invariant distribution volume of glucose in the accessible compartment.

The following equalities are used on the basis of the assumptions of tracer indistinguishability: (4) (5) i.e., the fractional clearance of the EGP-mimicking glucose tracer determines the fractional clearance of the EGP-originating glucose and, similarly, the fractional clearance of the meal absorption-mimicking glucose tracer determines the fractional disposal of meal-originating glucose. Under the assumption that the 2CMs correctly represent the glucose kinetics, *k*_{01,EM} = *k*_{01,MM}. The difference between estimates of these two fractional clearances indicates the extent of model misspecification.

The concentration of the total glucose in the plasma sample is determined by a direct measurement. The concentration of tracer species S = EM, M, and MM in the sample is determined as (6) where TTR_{S} represents the tracer-to-tracee ratio of glucose tracer *S* over the native glucose in the sample. TTR_{S} is determined by the gas chromatography-mass spectrometry (GC-MS) or liquid chromatography-mass spectrometry (LC-MS) using calibration curves or, in case of spectra overlap and tracer recycling, an isotopomer correction method (11, 15, 20). The ratio on the right-hand side of *Eq. 6* represents the concentration of native glucose originating from the EGP and the native glucose in the meal.

The concentration of the native glucose originating solely from the EGP is obtained with the use of a model-independent formula (7) where D_{N} and D_{M} represent the amount of native and labeled glucose in the meal, respectively.

The time-invariant model parameters are assumed to be identical across all glucose species with *k*_{21} = 0.05/min, *k*_{12} = 0.07/min, and V_{1} = 160 ml/kg (11). The forcing functions *u*_{EM} and *u*_{MM} are assumed to be known with absolute accuracy. The unknowns to be estimated are the EGP (= *u*_{EN}), the R_{a meal} (= *u*_{M}), and the R_{d} (= *k*_{01,T} Q_{1,T}).

#### Traditional solution: model transformation.

Following the work by Steele (16), the traditional solution (1) utilizes data-determined time derivatives of glucose concentration and concentration ratios while transforming *Eqs. 1*–*3* to calculate directly the required quantities.

The time variant *k*_{01,EM} associated with the EGP-mimicking tracer can be derived from *Eqs. 1* and *3* as (8) Substituting *k*_{01,EM} from *Eq. 8* into *Eq. 1* for S = E while utilizing the identity given by *Eq. 4* leads to (9) Similarly, *k*_{01,MM} can be obtained as (10) and the R_{a meal} obtained as (11) Recalling that the R_{d} is defined as *k*_{01,T} Q_{1,T} and that (12) the R_{d} can be calculated from *Eqs. 1* and *3* for S = T as (13) The calculation of the EGP according to *Eq. 9* employs the concentration ratio G_{EM}/G_{E}, its time derivative d(G_{EM}/G_{E})/d*t*, and the concentration G_{E}. The values of G_{EM}/G_{E} and G_{E} are obtained by using a smoothing approach such as the optimized optimal segments program (2) or the more theoretically sound stochastic regularization method (5). The smoothing approach interpolates data between measurements to facilitate the provision of equidistantly spaced concentration ratios at time points T_{i}, *i* = 1…N, with d*t* = T_{i +1} − T_{i} normally equal to 5 min (12). Experimental data are usually nonequidistantly spaced with more frequent sampling around the experimental perturbation (note that the experimental sampling schedule given by *t*_{k}, *k* = 1…M, is a subset of the calculation sampling schedule given by T_{i}). Irrespective of the smoothing method, the time derivative is assumed constant in the period T_{i} to T_{i +1} and is normally obtained as [(G_{EM}/G_{E})_{i+1} − (G_{EM}/G_{E})_{i}]/(T_{i +1} − T_{i}), where (G_{EM}/G_{E})_{i} is the smoothed concentration ratio at time T_{i}. Smoothing is essential, as calculations without smoothing lead to nonphysiological oscillations due to ill conditioning. The masses Q_{2,EM} and Q_{2,E} in the last expression on the right-hand side of *Eq. 9* are approximated using a recursive formula, which assumes a piecewise linear G_{EM} and G_{E} between T_{i} and T_{i +1}, facilitating an efficient implementation in a spreadsheet (12).

A similar approach is used to calculate the R_{a meal} from *Eq. 11*. In this case, the concentration ratio G_{MM}/G_{M}, its time derivative, and the concentration G_{M} are obtained by smoothing and interpolating smoothed data to evaluate the concentration ratios between measurements. The recursive calculations of Q_{2,MM} and Q_{M} are also required.

The calculation of the R_{d} using *Eq. 13* utilizes the concentration G_{T} and its time derivative. These are again smoothed and interpolated at regular time intervals. The recursive calculations of Q_{2,T} are needed.

Overall, the calculations of the EGP, the R_{a meal}, and R_{d} require the evaluation of the following indexes at regular time points: G_{EM}/G_{E}, G_{E}, G_{MM}/G_{M}, G_{M}, and G_{T}. This involves smoothing and interpolating between sampling times. Time derivatives of G_{EM}/G_{E}, G_{MM}/G_{M}, and G_{T} are required. The calculations approximate surrogates of masses Q_{2,EM}, Q_{2,E}, Q_{2,MM}, Q_{M}, and Q_{2,T}.

Although numerically efficient, the calculation procedure suffers from three shortcomings. First, smoothing of the concentration ratios G_{EM}/G_{M} and G_{MM}/G_{M} distorts the characteristics of the measurement error. This is of particular concern for G_{MM}/G_{M} following meal digestion because concentrations of both G_{MM} and G_{M} are small and the ratio is greatly influenced by the measurement error. Basu et al. (1) suggested overcoming this problem by omitting the first G_{MM}/G_{M} ratio obtained at 10 min from data processing.

Second, smoothing regularizes the concentrations and concentration ratios but not the underlying clearance rates *k*_{01,S} and metabolic fluxes such as the EGP the R_{a meal}, which, on physiological basis, should conform to the regularity assumption. Due to experimental constraints, infusion rates *u*_{EM} and *u*_{MM} are normally piecewise constant, and when a step change occurs the resulting plasma concentration of *G*_{EM} and *G*_{MM} demonstrate a local shape change, with the time derivative of *G*_{EM} and *G*_{MM} presenting a step change. However, these shape changes will be smoothed out by the smoothing algorithm when evaluating G_{MM}/G_{M} and G_{EM}/G_{E} and particularly their time derivatives, which enter *Eqs. 9* and *11*. These two equations also utilize the discontinuous infusion rates *u*_{EM} and *u*_{MM} and may return, in consequence, jagged EGP and R_{a meal}. This will also cause the R_{d} to be jagged as both the EGP and R_{a meal} enter *Eq. 13*. In summary, smoothing plasma concentrations and ratios of plasma concentrations when using discontinuous tracer infusions result in a jagged and thus nonphysiological EGP, R_{a meal}, and R_{d}.

Finally, numerical errors are introduced through an approximate solution of *Eqs. 1* and *2*. The recursive calculation of Q_{2,S} surrogate, which assumes a piecewise linear Q_{1,S} between T_{i} and T_{i +1}, whereas in fact Q_{1,S} is a of sum of two exponentials plus a linear component (see *Eq. A9* in the appendix), introduces a numerical error, which would be of particular importance during dynamic conditions. An additional numerical error is introduced by assuming a constant value of Q_{2,S} surrogate during the calculation interval. A short calculation interval (i.e., 5 min) can reduce the impact of these numerical errors.

*Eqs. 9* and *11* are excellent in demonstrating the benefit of using *u*_{EM} and *u*_{MM}, which mimic the EGP and R_{a meal}. Then the calculations become virtually model independent. However, due to intersubject variability, constant ratios G_{EM}/G_{M} and G_{MM}/G_{M} are difficult to achieve on individual basis, and other methods to calculate the EGP, R_{a meal}, and R_{d} may be more appropriate.

#### New calculation method: model retention combined with maximum likelihood and regularization.

The method utilizes the maximum likelihood (ML) theory (7) combined with the regularization of Tikhonov et al. (18) to achieve a theoretically coherent computational framework. Under a mild assumption of the constancy of *k*_{01,S} between the measurement points, this approach also avails itself for implementation in a spreadsheet without the need for numerical approximations.

Briefly, the method estimates the time-variant *k*_{01,EM} using the plasma measurements of the EGP-mimicking tracer. Then, the *k*_{01,EM} estimate is used to calculate the EGP from the plasma measurements of the endogenous glucose. A similar approach is used to estimate the R_{a meal}. First, *k*_{01,MM} is estimated from plasma measurements of the meal absorption-mimicking tracer. Then, the R_{a meal} is estimated from the plasma measurements of the tracer included with the meal. Finally, *k*_{01,T} is estimated from the total plasma glucose concentration, utilizing the EGP and R_{a meal} estimates and yielding the R_{d}.

The detail derivation follows. Let Ḡ_{EM,k}, *k* = 1*…M*, represent the measurement of the EGP-mimicking tracer at time *t*_{k}, *t*_{k +1} > *t*_{k}, with the measurement error normally distributed, uncorrelated, with a zero mean and a standard deviation σ_{EM,k}. The logarithm of the likelihood function LLF(EM) defining the objective function for the estimation of *k*_{01,EM}(*t*) consists of two components: (14) The first component, LLF_{fit}(EM), evaluates the model fit (15) where G_{EM} (*t*_{k}) is the solution of *Eqs. 1*–*3* for a piecewise constant *k*_{01,EM}(*t*), defined as (16) with *k*_{01,EM,k} > 0, *k* = 1… M − 1.

The second component, LLF_{reg}(EM), evaluates regularity of *k*_{EM}(*t*). The variable is assumed to be a Wiener process, where ω_{EM} represents the weighting factor balancing the relative importance of the LLF_{reg} against the LLF_{fit} component. The function ξ(*t*) allows smoothness of *k*_{01,EM}(*t*) to be reduced at the time of an experimental perturbation such as the meal intake, utilizing the knowledge of the experimental design in the calculation process. Assuming piecewise constant fractional clearance *k*_{01,EM}(*t*), the differences *k*_{01,EM,k +1} − *k*_{01,EM,k} are thus normally distributed, uncorrelated, with zero mean and a standard deviation ξ_{k}ω_{EM} (17) The adoption of the Weiner process expresses our belief that changes in *k*_{01,EM}(*t*) over a short time period are more likely than proportional changes over a longer time period, allowing variations in *k*_{01,EM}(*t*) to be captured when frequent sampling is used.

The variables ξ_{k} are defined as (18) The constant ξ,ξ > 1, is suitably chosen. In the present study, ξ = was used.

The LLF_{reg}(EM) is then defined as (19) Maximizing the likelihood is identical to minimizing its negative value, and thus the solution is found as (20) where Q_{20,EM} in *Eq. 2* is obtained assuming steady-state conditions: (21) Differentiating LLF(EM) with respect to ω_{EM} in *Eq. 20*, setting the derivative to zero, and solving for ω_{EM} gives the solution at which the minimum of *Eq. 20* is attained (22) Substituting for ω from *Eq. 22* into *Eq. 19* and removing constants, our minimization problem becomes (23) The weighting factor ω_{EM} is absent in *Eq. 23*. The regularity of *k*_{01,EM} is fully determined by the standard deviations of the measurement error σ_{EM,k}. A high-measurement error will enable only overly smooth *k*_{01,EM} to be determined, and a low-measurement error facilitates a refined estimation of *k*_{01,EM}.

The estimation of the EGP follows a similar path. The difference is that the EGP is continuous piecewise linear defined by a sequence, *u*_{E,k}, *k* = 1…M: (24) The minimization problem defining solution for the EGP is then (25) where G_{E}(*t*_{k}) is the solution of *Eqs. 1*–*3* for *k*_{E}(*t*) obtained from *Eq. 23*. Ḡ_{E,k} represents the *k*th measurement of the endogenous glucose and σ_{E,k} the standard deviation of the associated measurement error. Q_{20,E} is again obtained assuming steady-state conditions: (26) The setup of the problem to estimate *k*_{01,MM}(*t*) and *u*_{M}(*t*) parallels that of *k*_{01,EM}(*t*) and *u*_{E}(*t*). In *Eqs. 15*–*25*, terms with subscripted EM are replaced by terms with subscripted MM, and terms with subscripted E are replaced by terms with subscripted M. The only exception is that the initial conditions Q_{10,MM} and Q_{10,M} are dropped from the optimization problems, as they are assigned zero value given that the appearance of the meal tracer and the meal absorption-mimicking tracer commence after the time origin. The R_{a meal} is obtained as (27) Finally, the estimation of the R_{d} is set up as the estimation of *k*_{01,T} for consistency with the estimation of *k*_{01,EM} and *k*_{01,MM} and reflecting physiological considerations that the fractional clearance rather than the glucose flux is regular. The optimization problem is defined by *Eqs. 14*–*23* by substituting terms with subscripted EM by terms with subscripted T. The initial condition Q_{10,T} is dropped from the optimization problem defined by *Eq. 23* since, for consistency, Q_{10,T} is evaluated from the initial conditions for the endogenous glucose and the EGP-mimicking tracer: (28) The appearance rate *u*_{T}(*t*) is defined by *Eq. 12*. Given optimized *k*_{01,T}(*t*), the R_{d} is obtained as (29) The minimization problems in *Eqs. 23* and *25* require solution of the 2CM defined by *Eqs. 1*–*3* for a piecewise constant *k*_{01,S}(*t*) and a piecewise linear *u*_{S}(*t*). A closed-form solution exists (see appendix).

The estimation of the EGP, R_{a meal}, and R_{d} is set up as a sequence of five optimization problems. The five optimization problems estimate *1*) three nonnegative sequences, *k*_{01,EM}, *k*_{01,MM}, and *k*_{01,T}, *k* = 1…M − 1; *2*) two sequences, *u*_{E} and *u*_{M}, *k* = 1…M; and *3*) two initial conditions, Q_{10,EM} and Q_{10,E}. The measurements Ḡ_{EM,k}, Ḡ_{E,k}, Ḡ_{MM,k}, Ḡ_{M,k}, Ḡ_{T,k}, *k* = 1…M, and associated standard deviations of the measurement error, σ_{EM,k}, σ_{E,k}, σ_{MM,k}, σ_{M,k}, and σ_{T,k}, are utilized in the estimation process. The results of the optimization problem for *k*_{01,EM} feed into the optimization problem for the EGP. Similarly, the results of the optimization problem for *k*_{01,MM} feed into the optimization problem for the R_{a meal}. Finally, the optimization problem to estimate *k*_{01,T} utilizes the EGP and R_{a meal}.

These five optimization problems can be replaced by one large optimization problem, summing up the five objective functions defined for the individual subproblems; i.e., the overall likelihood function LLF is defined as (30) and the optimization problem consists of minimizing the negative of the LLF: (31) This facilitates a more comprehensive utilization of the experimental data with measurements of the endogenous glucose feeding into the estimation of *k*_{01,EM}, for example. The drawback is an increased computational complexity. In the present work, the five optimization problems were first solved sequentially, and then *Eq. 31* was used to obtain the final estimate of the unknown variables.

The calculations were implemented in a spreadsheet utilizing the “solve” function to solve the nonlinear minimization problems. A sample spreadsheet demonstrating the calculations can be obtained from the corresponding author free of charge by academic institutions and adopted for noncommercial projects.

#### Sample experimental study.

A 30-yr-old healthy female (weight 62 kg, body mass index 25.7 kg/m^{2}) was studied after overnight fast. The subject received a primed, constant, continuous infusion of 45.5 nmol·kg^{−1}·min^{−1} [1,2,3,4,5,6,6-^{2}H_{7}]glucose (the EGP-mimicking glucose tracer; Cambridge Isotope Labs) for 120 min prior to the digestion of a liquid mixed meal (600 kcal, 75 g of glucose, 32.9 g of protein, and 20 g of fat) containing 3.11 g of [2-^{2}H_{1}]glucose (Omacron). Following meal digestion, the infusion of [1,2,3,4,5,6,6-^{2}H_{7}]glucose was piecewise linear in a fashion anticipating the time-varying profile of the EGP. A piecewise linear infusion of [U-^{13}C;1,2,3,4,5,6,6-^{2}H_{7}]glucose (the meal absorption-mimicking glucose tracer; Cambridge Isotope Labs) was commenced 10 min after meal digestion and reached a peak of 5.43 nmol·kg^{−1}·min^{−1}. Figure 1 shows the time profile of the EGP-mimicking glucose tracer and the meal absorption-mimicking glucose tracer.

Samples were taken at −30, −15, −5, 0, 10, 20, 30, 40, 50, 60, 70, 90, 120, 150, 180, 210, 240, 270, 300, 330, and 360 min, analyzed for glucose, and subjected to GC-MS/LC-MS analysis to determine TTR_{s} of the three glucose tracers. In brief, plasma samples for GC-MS were derivatized after protein precipitation to the aldonitrile pentacetate with hydroxylamine hydrochloride-acetic anhydride. GC/electron impact-mass spectrometry analysis was performed on an Agilent model 6890/5973 with a 7673 autosampler (Agilent). Fragments representing carbons 1–5 or 5,6 were analyzed for [2-^{2}H_{1}]glucose and [1,2,3,4,5,6,6-^{2}H_{7}]glucose, respectively. [U-^{13}C;1,2,3,4,5,6,6-^{2}H_{7}]glucose was measured by LC-MS as previously described (14). Plasma glucose was measured using an enzymatic method Beckman Glucose Analyzer II (Beckman, Fullerton, CA). TTRs were obtained from isotope ratios using calibration curves, as isotopomer effect is negligible with the adopted tracers. Tracer recycling was below the limit of detection in constant infusion control experiments covering the relative ratios of plasma tracer concentrations and was therefore considered negligible.

#### Data analysis: traditional solution.

The traditional analysis adopted smoothing and *Eqs. 9*, *11*, and *13* to calculate the EGP, the R_{a meal}, and the R_{d}. Data processing involved the use of the CODE program, implementing stochastic regularization technique (9) to obtain G_{EM}/G_{E}, G_{E}, G_{MM}/G_{M}, G_{M}, and G_{T} at 5-min intervals. A constant coefficient of variation for each quantity was adopted to determine the extent of smoothing. The ratio G_{MM}/G_{M} at 10 min was excluded from the analysis as suggested previously (1). Time derivatives were also calculated at 5-min intervals as the difference between interpolated values divided by the 5 min. The recursive technique proposed by Mari (12) was used to obtain the surrogate values of Q_{2,EM}, Q_{2,E}, Q_{2,MM}, Q_{M}, and Q_{2,T}.

#### Data analysis: new method based on maximum likelihood.

The five optimization problems defined by *Eqs. 23* and *25* were, in the first instance, solved sequentially. Then, *Eq. 31* was used to obtain the final estimate of the EGP, R_{a meal}, R_{d}, *k*_{01,EM}, *k*_{01,MM}, and *k*_{01,T}. The measurement error associated with the total glucose was assumed to be multiplicative with 1.5% CV. The SD of the measurement error associated with the endogenous glucose was assigned values identical to that associated with the total glucose, i.e., σ_{E,k} = σ_{T,k}, reflecting the calculation methods of the endogenous glucose. The standard deviations of the measurement error associated with the TTR_{s} of [1,2,3,4,5,6,6-^{2}H_{7}]glucose, [U-^{13}C;1,2,3,4,5,6,6-^{2}H_{7}]glucose, and [2-^{2}H_{1}]glucose were set to 1.32 × 10^{−5}, 2.40 × 10^{−6}, and 9.68 × 10^{−4} (unitless), respectively. The standard deviations σ_{EM,k}, σ_{MM,k}, and σ_{M,k} were calculated from σ_{T,k} and the standard deviations associated with the respective TTR using the error propagation technique, adopting *Eq. 6*. The optimization problem defined by *Eq. 31* was implemented in a spreadsheet.

#### Evaluation using simulation study.

A comparison of the two approaches employed synthetic data sets generated by a glucoregulatory model. The synthetic profiles were analyzed with the traditional and new methods to reconstruct EGP, R_{a meal}, and R_{d} profiles. The differences between the actual and reconstructed EGP, R_{a meal}, and R_{d} profiles were summarized using the root mean square error (RMSE).

Two categories of three-tracer synthetic experiments were generated with six experiments in each category. In the first category, the time-varying infusions of the EM and MM tracers were identical in all six experiments and were based on the information from the literature (1). In the second category, the time-varying infusions of the EM and MM tracers were optimized using individual EGP and R_{a meal} profiles. This allowed errors because of the experimental nature and due to the model misspecification to be separated.

#### Generation of synthetic data sets.

The synthetic data sets were generated with a glucoregulatory model consisting of 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 used a two-compartment structure, as described by Hovorka et al. (8). The submodel of insulin secretion assumed a linear relationship between plasma glucose and insulin secretion (10). A one-compartment model of the insulin kinetics was adopted (8). A submodel of the glucose kinetics assumed a two-compartment structure with three insulin action compartments representing insulin action on the glucose distribution/transport, disposal, and EGP suppression (8).

Six parameter sets were randomly generated from prior distributions (8, 10) to represent six synthetic subjects. Following the ingestion of 75-g carbohydrate meal, these parameter sets resulted in the time to peak of the gut absorption of 68 ± 9 min (mean ± SD) and the maximum suppression of the EGP of 81 ± 6% at 84 ± 30 min postmeal.

Two categories of three-tracer simulated experiments were executed. In the first category, the EM and MM tracers were given as piecewise constant infusions, using shapes defined by Basu et al. (1). In the second category, the EM and MM tracers were given again as piecewise constant infusions, with the time resolution defined by Basu et al., but the profiles were obtained by sampling the individual EGP and the individual gut absorption at midpoints of the time resolution scheme. In principle, the second-category experiments should make the estimation of the EGP, R_{a meal}, and R_{d} more accurate as the EM and MM tracer infusions more closely follow the EGP and the gut absorption, except for the piecewise approximation of the EGP and gut absorption profiles.

The model-derived total plasma glucose concentration and the TTR of the EM, MM, and M tracers were sampled with the sampling schedule described in *Sample experimental study*. A measurement error at the extent described in *Sample experimental study* was added to the measurements.

#### Data analysis.

The noisy measurements were subjected to the analysis by the traditional and the new method, as described in *Sample experimental study*. This provided estimates of the EGP, R_{a meal}, and R_{d}.

#### Error assessment.

The difference between the actual and the estimated EGP, R_{a meal}, and R_{d} profiles was summarized by evaluating the RMSE at a 10-min interval, i.e. (32) where r_{a,i} and r_{e,i} represent the actual and estimated profiles, respectively, sampled at time *t*_{i} = (*i* − 1) × 10 min.

#### Statistical analysis.

A two-way analysis of variance assessed the difference in the RMSE of the EGP, R_{a meal}, and R_{d} with the method (new vs. traditional) being one factor and the experimental category (literature-based vs. individually optimized EM and MM tracer infusions) being the other factor. A paired *t*-test was used to evaluate the difference between the RMSE associated with the R_{a meal} calculated with the traditional method comparing the effect of the literature-based vs. individually optimized EM and MM infusions. A similar, paired *t*-test was used to evaluate the R_{a meal} obtained by the new method, the R_{d} obtained by the traditional method, and the R_{d} obtained by the new method, as an interaction existed between the method and the experimental category when evaluating the RMSE associated with the R_{a meal} and R_{d}.

## RESULTS

#### Sample experimental study.

The traditional analysis adopting smoothing and *Eqs. 9*, *11*, and *13* gave the EGP, the R_{a meal}, and the R_{d}, as shown in Fig. 2. The corresponding fluxes obtained with the use of the ML method are also shown in Fig. 2.

The sample study indicates a nearly identical EGP obtained with the two computational approaches. The R_{a meal} and the R_{d} demonstrate considerably higher swings in the late postprandial period with the traditional method. Fast but low-amplitude oscillations in R_{d} are also observed during the early postprandial period with the traditional method.

The concentrations of the three tracers are shown in Fig. 3. The plots also show the reconstructed tracer concentrations using the two computational approaches. The fit to the EGP-mimicking tracer is nearly identical, reflecting a nearly identical EGP estimate by the two methods. Differences in the fit to the data between the two approaches are observed with the meal-mimicking tracer and the tracer included with the meal in the second part of the study.

The endogenous glucose and plasma (total) glucose concentration together with reconstructed profiles are shown in Fig. 4.

The concentration ratios G_{EM}/G_{M} and G_{MM}/G_{M}, which play the dominant role with the traditional approach, are shown in Fig. 5. The reconstructed ratios are also shown. The anticipated EGP and R_{a meal} differed substantially from the calculated rates (see Figs. 1 and 2), highlighting the difficulties in achieving constant concentration ratios on an individual basis.

The fractional clearance rates *k*_{01,EM}, *k*_{01,MM}, and *k*_{01,T} are shown in Fig. 6. The large value of *k*_{01,MM} calculated by the traditional approach at 20 min is caused by a nonzero value of G_{M} at 20 min. The differences between *k*_{01,MM} by the two methods in the second part of the study are linked to the differences in fitting the measured G_{MM} concentration; see Fig. 3, *middle*. The two plots indicate the extent of ill conditioning. Relatively large differences in *k*_{01,MM} result in small differences in G_{MM}.

#### Evaluation using simulation study.

RMSEs associated with the new and the traditional methods and two experimental categories are shown in Table 1.

The results document a significantly better performance of the new method when calculating the EGP, R_{a meal}, and R_{d}. The EGP was estimated with the lowest RMSE, followed by R_{a meal} and R_{d} by both methods.

The optimization of the EM and MM tracer infusions reduced the RMSE associated with the EGP by both methods. The optimization of the tracer infusions also reduced the RMSE associated with the R_{a meal} and R_{d} using the new method (*P* < 0.05, paired *t*-test). Unexpectedly, the optimization resulted in an increase of the RMSE associated with the R_{a meal} and R_{d} using the traditional method (*P* < 0.001 and *P* < 0.05, paired *t*-test). On a visual inspection (data not shown), the increase in the RMSE was attributable to a marked mismatch between the actual and estimated R_{a meal} in the first 30 min. The traditional method overestimated the actual R_{a meal} with an accelerated time to peak of the gut absorption at 20 min, whereas the true time to peak was 70 min. Further overestimation to a lesser extent and a jagged R_{a meal} profile followed from 60 min onward. The overestimation of the R_{a meal} propagated into the overestimation of the R_{d}, as the mass balance had to be maintained.

## DISCUSSION

A new method is proposed to calculate the EGP, R_{a meal}, and R_{d} during a three-tracer study design. The method retains the two-compartment formulation of the glucose kinetics and avoids the transformation of the measurement error. Physiologically justified assumptions are made about the smoothness of the underlying fractional clearance rates and glucose fluxes, supporting plausibility of the calculations.

The traditional approach divides the calculation process into two stages. First, data are smoothed and interpolated. Second, the smoothed data, assumed to be error-free, enter recursive formulae. The computational complexity resides within the smoothing stage, which includes a nonlinear optimization problem (2, 5). The new approach contains only one stage. The computational complexity is comparable to the smoothing stage of the traditional approach. The nonlinear optimization problem is integrated by placing smoothness assumption on the underlying metabolic indexes. Yet the computations can be implemented in a spreadsheet.

The new method has further benefits when used during a piecewise constant infusion of the EGP-mimicking and meal-mimicking tracers. Step changes in the infusion rate result in “jumps” in time derivative of the tracer concentration, but the jumps are filtered out by the smoothing stage of the traditional approach. However, as the piecewise constant infusions also enter the recursive formulae, a jump (discontinuity) is introduced in the calculated EGP and R_{a meal}. The R_{a meal} is particularly affected, as the meal-mimicking tracer infusion varies extensively. Optimally, the discontinuity in the infusion rate and the measured concentration should both enter the recursive formulae and cancel out each other.

With the traditional approach, the calculation of the R_{d} and R_{a meal} is sequential. The estimated R_{a meal} enters the formula to calculate the R_{d}; see *Eqs. 12* and *13*. The drawback is that any error in the R_{a meal} will be accompanied by a proportional error in the R_{d}. This is apparent in Fig. 2, *middle* and *bottom* [Mari method (12)]. The oscillations in the R_{a meal} from 210 min until the end of the experiment are followed by oscillations in the R_{d} of a similar magnitude.

The propagation of the estimation error is a consequence of the sequential nature of the calculations. Errors in both the EGP and R_{a meal} will introduce a comparable error in the R_{d}, although the extent of the former is expected to be smaller given the lower anticipated extent of the error associated with the EGP.

The new approach eliminates the propagation of the error. The calculation is no longer sequential, and all fluxes and fractional clearance rates are calculated in one go; see *Eq. 31*. Conceptually, this can be visualized as information sharing. All five measurements, the three tracer concentrations, the endogenous glucose concentration, and the total glucose concentration influence the calculation of the EGP, R_{a meal}, and R_{d}. The assumption of smoothness of the R_{d} forces the R_{a meal} to reduce its oscillatory pattern (see Fig. 2, *middle*), ML method, albeit at the expense of a small misfit to the [U-^{13}C;1,2,3,4,5,6,6-^{2}H_{7}]glucose (Fig. 3, *middle*) but an improved fit to the [2-^{2}H_{1}]glucose (Fig. 3, *bottom*). Although less oscillatory, the R_{d} provides a better fit to the total glucose; see Fig. 4, *bottom*.

Another drawback of the traditional approach is the transformation of the measurement error. The first 20–30 min following meal ingestion provide low values of G_{MM} and G_{M}. The measurement error will represent a considerable proportion of the two measurements and will impact on the ratio G_{MM}/G_{M}. The traditional approach solves this problem by excluding the ratios G_{MM}/G_{M} during the first 10 min from smoothing. However, in principle, this problem perseveres, although to a lesser extent, for other G_{MM}/G_{M} ratios and also G_{EM}/G_{E} ratios. In our particular example, the effect is minimized by the adoption of highly precise GC-MS and LC-MS measurement techniques. The new computational approach avoids using the ratios and utilizes solely measured concentrations for which the measurement error is determined on a coherent basis.

The fractional turnover rate *k*_{01} differs when estimated using the EM and MM tracers; see Fig. 6. The differences in *k*_{01} estimates reflect the extent of model misspecification and, presumably, to a smaller extent, the effect of the measurement error on parameter estimation. Conceptually, if a model is incorrect, different inputs, such as in our case different EM and MM tracer infusions, will provide different parameter estimates. If the model is a faithful representation of the reality, the parameter estimates that are not influenced by the input bar the effect of the measurement error.

In principle, glucose cycling and nonequivalent tracer loss could also contribute to the observed differences in *k*_{01}. However, our assessment of glucose cycling indicates that this is unlikely. The [1,2,3,4,5,6,6-^{2}H_{7}]glucose was monitored as the C_{1}–C_{5} fragment containing four deuteriums (C_{2}–C_{5}). The enrichment of the [1,2,3,4,5,6,6-^{2}H_{7}]glucose molecule was five- to 10-fold less (1:5–10) than the enrichment of the [2-^{2}H_{1}]glucose for most of the experiment, being close to 1:3 only at the very beginning and the extreme end of the meal. As expected, in test infusions we found no evidence of tracer recycling of the [1,2,3,4,5,6,6-^{2}H_{7}]glucose onto [2-^{2}H_{1}]glucose because the fragment would have added an M + 2 or M + 3. Similarly, with the [U-^{13}C;1,2,3,4,5,6,6-^{2}H_{7}]glucose, we found no measurable evidence of recycling onto the [1,2,3,4,5,6,6-^{2}H_{7}]glucose in test infusions at the expected ratios during the study.

Although less transparent, the new method also benefits from the constant enrichments, i.e., a zero time derivative of the respective ratios. To illustrate this point, let us consider the calculation of the EGP by the Mari approach (12). *Eq. 9* includes the time derivative of the ratio of the tracer to endogenous glucose, and the minimization of the derivative is achieved by mimicking the EGP with the EM tracer. The derivation of the Mari formula to calculate the EGP involves the substitution of the formula calculating *k*_{01}, *Eq. 8*, into *Eq. 1* and associated algebraic manipulations. The point here is that *k*_{01} is an “intermediate” product, parts of which conveniently drop out in the Mari model. The benefit is a simplification of the formula and a clear exposition of the effect of minimizing the time derivative. The drawback is that physiological information about the smoothness of *k*_{01} is lost and is replaced by assuming smoothness of the ratio MM tracer to the endogenous glucose.

In the new method, the intermediate step to calculate *k*_{01} is retained. This may confound the understating that the minimization of the time derivate also reduces the effect of model misspecification, but this principle holds here, too. Unlike algebraic manipulations adopted by the Mari model, the new method uses a numerical estimate of *k*_{01}. The benefit is the use of the prior information on *k*_{01} smoothness. It therefore follows that, if the infusions of the EM and MM tracers exactly mimic the EGP and the gut absorption, the new method will provide accurate estimates of the EGP and the gut absorption despite the different *k*_{01} values estimated from the EM and MM tracers.

The Mari method and, similarly, the new method assume that the fractional transfer rates *k*_{21} and *k*_{12} are time variant and set to physiologically feasible values and identical to those adopted by Basu et al. (1). However, it is known that these fractional rates are stimulated by insulin with a greater effect of insulin on *k*_{21} than on *k*_{12} (4), pinpointing one source of model misspecification that will contribute, when the EM and MM tracers differ from the EGP and the gut absorption, to the computational errors in the EGP, R_{a meal}, and R_{d}. However, these errors will be smaller than those incurred when using a single compartment model (1).

The volume of distribution V_{1} was set at value of 160 ml/kg, comparable to the 146 ml/kg used by Mari (12) and slightly higher than 130 ml/kg used by Basu et al (1). It is identical to the value obtained during modeling of the intravenous glucose tolerance test with a two-compartment model (11).

The theoretical expectations of superiority of the new method were confirmed in a simulation study. The new method outperformed significantly the traditional method when estimating the EGP, R_{a meal}, and R_{d} with the greatest RMSE improvement associated with the calculation of R_{a meal}. Furthermore, the new method benefited from the optimized infusion profiles of the EM and MM tracers. The optimized infusions resulted in halving the RMSE associated with R_{a meal}. Smaller relative improvements were observed in the RMSE associated with the R_{d} and EGP.

Unexpectedly, the optimized tracer infusions resulted in a deteriorated RMSE associated with the R_{a meal} and subsequently the R_{d} using the traditional method. This is attributable to a marked overestimation of the R_{a meal} in the first 30 min, although additional overestimation followed from 60 to 120 min with jagged profile throughout. This jagged profile was not observed with the EGP, suggesting that the calculations of the R_{a meal} and R_{d}, but not the EGP, suffer from the piecewise constant infusion of the tracers with the traditional method.

The reasons for the overestimation of the R_{a meal} in the first 30 min may originate in the effect of the measurement errors on the calculations. With the optimized MM tracer infusion, the infusion rate of the MM tracer is considerably smaller in the initial part of the experiment. This smaller infusion rate will result in a lower TTR of the MM tracer. Given that the measurement error associated with the TTR of the MM tracer is additive, the signal-to-noise ratio at the early part of the experiment will be lower with the optimized infusion. These highly noisy measurements are used as a numerator when evaluating and smoothing the G_{M}/G_{MM} ratio. It appears that the smoothing process of highly noisy measurements is the culprit of the problem.

The computations associated with the new approach can be executed in a spreadsheet, given the explicit solution of the two-compartment model; see appendix. The solution utilizes an assumption that the fractional clearance rate *k*_{01,S} is constant between time instances. Figure 6 shows the time evolution of the fractional clearance rates, indicating jumps in the fractional clearance rate associated with *k*_{01,MM}. The jumps could be reduced by using a finer time resolution.

The new computational approach can be used in different experimental scenarios. This includes the two tracer study designs to calculate the EGP, R_{a meal}, and R_{d} or the glucose clamp technique to calculate the EGP and R_{d}. *Eq. 31* needs to be suitably modified to represent the reduced experimental complexity. All other components remain unchanged.

In conclusion, a new computational approach has been developed, increasing feasibility of the EGP, R_{a meal}, and R_{d} calculated from data collected during the three-tracer experiment. The new approach uses more physiologically justified assumptions and treats coherently the measurement error.

## APPENDIX

Here we describe the closed-form solution of the glucose kinetics model defined by *Eqs. 1* and *2*. The glucose species subscript *S* is dropped from the notation, as the solution is identical for all species.

The solution is recursive. It assumes that *k*_{21} and *k*_{12} are time invariant. Given a possibly nonequidistant time sequence *t*_{k}, *k*= 1…M, it is assumed that *k*_{01}(*t*) is nonnegative, piecewise constant, i.e., *k*_{01}(*t*) = *k*_{01,k} for *t*_{k} *t*_{k+1}, *k*_{01,k} > 0. It is further assumed that *u*(*t*) is continuous, piecewise linear, *u*(*t*) = *u*_{k} + (*u*_{k+1} − *u*_{k})/(*t*_{k+1} − *t*_{k})(*t*−*t*_{k}) for *t*_{k} ≤ *t* < *t*_{k+1}.

Given the amounts Q_{1}(*t*_{k}) and Q_{2}(*t*_{k}), it can be shown that amounts Q_{1}(*t*_{k+1}) and Q_{2}(*t*_{k+1}) can be obtained as (A1) (A2) (A3) (A4) (A5) (A6) (A7) (A8) (A9) (A10) where r(*t*_{k}) is the rate of change in *u*(*t*_{k}) and A(*t*_{k}), B(*t*_{k}), C(*t*_{k}), D(*t*_{k}), E(*t*_{k}), F1(*t*_{k}), and F2(*t*_{k}) are auxiliary variables.

## GRANTS

Support by the European Foundation for the Study of Diabetes, the European Union-funded Clincip Project (IST-2002-506965), National Institute of Diabetes and Digestive and Kidney Diseases Grant RO1-DK-61641, National Center for Research Resources Grant MO1-RR-12248, and the American Diabetes Association is gratefully acknowledged.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2007 by American Physiological Society