ogttMetrics: Data structures and algorithms for oral glucose tolerance tests

Informal derivation of insulin sensitivity via the minimal model

Following Bergman et al. (1979), let G(t) and I(t) denote time-dependent plasma concentrations of glucose and insulin respectively. Various time-dependent factors affect the trajectories of these concentration functions, and we will assume that derivatives of these function with respect to time, and partial derivatives of these functions with respect to relevant time-dependent variables, can be defined. Let Ġ denote the rate of change of glucose concentration in plasma over time. Glucose effectiveness is defined as E = –∂Ġ/∂G. This is described as “the quantitative enhancement of glucose disappearance due to an increase in the plasma glucose concentration” (Bergman et al., 1979, p. E673). At steady state, insulin sensitivity is S_I = ∂E/∂I. A four compartment model (model VI of Bergman et al. (1979)) leads to differential equations

                                                G′(t) = (p₁ – X(t))G(t) + B₀

and

                                                X′(t) = p₂X(t) + p₃I(t)

where X(t) is an abstract time-dependent function representing insulin action, G(t) is the time-dependent function representing glucose concentration, I(t) represents time-dependent insulin concentration, and B₀ represents “glucose balance” (difference between rates of hepatic release to circulation and uptake in peripheral tissue) extrapolated to zero glucose concentration. By the definition of glucose effectiveness, the first differential equation implies E(t) = X(t) – p₁, and, at steady state, X_SS = –I_SSp₃/p₂. This final expression is substituted into the expression for E just obtained, and after formal partial differentiation by I, we obtain S_I = –p₃/p₂.

A formal specification

In the following,

We consider the specific formalism for the dalla Man et al. minimal model given by Burattini et al. (2006):

with initial conditions G(0) = G_b and X(0) = 0.

Estimation

The procedure of Dalla Man et al. (2002) involves two phases. In the first phase, the system of ordinary differential equations (ODE) above is solved on the basis of provisional settings of unknown parameters. The solution yields pointwise predictions of glucose concentrations Ĝ_t with t ranging over the sampling time course of the OGTT. In the second phase, parameters of the ODE system are updated using non-linear least squares. The phases are iterated until the sum of squared discrepancies ∑_t(G_t – Ĝ_t)² converges to a minimum. Inputs to the algorithm are measured time series of glucose and insulin concentrations, and individual body weight; other quantities, such as glucose effectiveness (S_G) fraction of ingested dose absorbed (FA), and volume of distribution (V) are taken as fixed constants, with values derived from results of other experiments.

Programming considerations

The kernel of fitOneMinMod in the ogttMetrics package is

    model <- function(t, Y, parameters) {
        with(as.list(parameters), {
            dy1 = -(Sg + Y[2]) * Y[1] + Sg * Gb + (ra(a1, a2,
                a3, t, BW, D, FA, DC)/V)
            dy2 = -p2 * Y[2] + p2 * SI * (Insulin(t) - Ib)
            list(c(dy1, dy2))
        })
    }

Here the interface to lsoda in the deSolve package is employed (Hindmarsh, 1983; Petzold, 1983). The formal variables Y[1] and Y[2] represent G(t) and X(t) respectively; ra() and Insulin() are specially defined functions that return, for any given time in the course of the OGTT, the rate of glucose appearance, and insulin concentration, respectively. The lsoda solver is invoked in the function mmsolfn, whose inputs a1, a2, a3 are free parameters of a piecewise linear model for Ra(t), the rate of appearance of glucose; input SI is the target quantity of interest, the measure of insulin sensitivity. The values of free parameters are obtained by minimizing the sum of squared differences between observed glucose g and values predicted by the ODE system for current values of the unknown parameters:

    mmsolfn = function(a1, a2, a3, SI) lsoda(c(Gb, 0), t, model,
        c(a1 = a1, a2 = a2, a3 = a3, SI = SI))
    fit = nls(g ~ mmsolfn(a1, a2, a3, SI)[, 2], start = nlsinit,
        trace = nlstrace, control = fullNLScontrol)

Additional quantities BW, D, FA, DC are used to implement the constraint of Dalla Man et al. (2002)

in which BW is participant body weight, D is the dose of glucose ingested, and FA is the fraction of ingested glucose that is actually absorbed; DC is constant that determines the rate of exponential decay of glucose concentration in plasma past minute 120.

For concreteness, Figure 1 displays all components of a minimal model fitted to a single 120 minute OGTT.

Data management and reporting

In practice, OGTT series can be collected according to different protocols and may include additional biomarkers such as c-peptide concentrations. For flexible data management and analysis, we adopted the data structure of the MultiAssayExperiment package of Bioconductor. We extended this structure in a class called ogttCohort, which includes metadata about timing of concentration measures. Each biomarker series for each individual is stored as a column of an R matrix, with rows and columns coordinated across assays. Arbitrary additional sample-level information can be linked to assay data. High level functions getMinmodSIs and addMatsuda120 fit the minimal model or compute the Matsuda index for each series, and append results to the data container. Because the minimal model may be time-consuming to fit, support is provided for parallel computation of multiple models. Use of a compact formal representation of all the OGTT data collected on a cohort simplifies creation of generic reports and visualizations. Figure 2 is based on the QCplots function, that can be applied to any ogttCohort instance. The top two panels display aspects of marginal (time-specific) distributions using boxplots. The bottom two panels are views of joint distributions of features and samples using the biplot methodology of Gabriel (1971). Calibrated outlier detection, proceeding under the assumption that the OGTT series are multivariate normal with a common mean vector and unspecified covariance matrix, can be conducted for glucose and insulin series separately, using mvOutliers. The procedure of Caroni & Prescott (1992) is used.

Figure 2. Output of QCplots for a sample of 50 observations from the OMNICarb study in the obaSamp object distributed with ogttMetrics.

Top two panels are time-specific boxplots, bottom two are biplots based on principal components analysis of the 50x2 7-dimensional vectors of glucose and insulin concentrations in the dataset.

Numerical considerations

Theodorakis et al. (2017) mentions that the standard (proprietary, closed source) software tool SAAM-II (Barrett et al. (1998)) failed to produce accceptable estimates of insulin sensitivity in over one-third of 106 samples. Similar difficulties were encountered in the OMNICarb study. These challenges motivated us to create an open source solution that would foster investigation of aspects of glucose and insulin series for which the minimal model fails to converge, and allow comparison of alternative metrics of carbohydrate metabolism on large datasets. Figure 3 displays the SIexplorer interactive interface. Given a collection of OGTT results in an ogttCohort structure, the SI vs Matsuda panel shows the association between estimated SI, Matsuda’s index, and convergence status of the Burattini et al. formulation of the minimal model. The display is made with transformed axes (log10 for SI, square root for Matsuda’s index). Negative estimates of SI are Winsorized to the smallest positive estimate observed in the data. Positive correlation between the indices is apparent, and the general trend appears to be obeyed for the majority of estimates of SI for which the dalla Man et al., algorithm does not converge.

Figure 3. SIexplorer interactive display of association of insulin sensitivity, Matsuda index, and minimal model convergence with default settings.

Application to a cross-over trial

We created the ogttMetrics package to analyze data from the OMNICarb study (Sacks et al. (2014)). This study involved over 150 overweight individuals (BMI > 25kg=m²) whose systolic blood pressure was in the interval 120–159 mmHg, or diastolic blood pressure in the interval 70–90 mmHg. Individuals with diagnoses of diabetes, cardiovascular disease, or chronic kidney disease were excluded. Four experimental diets were designed to provide contrasting values of overall carbohydrate content and glycemic index of foods consumed. Carbohydrate and glycemic index each had two levels denoted C and c (G and g) respectively, leading to the set (CG, Cg, cG, cg) of experimental diets. Each patient received a randomly ordered sequence of diets from this set, consuming each assigned diet for five weeks, with a pause of two weeks between diets. At the end of each feeding period a 120-minute OGTT protocol was administered. As noted previously, attempts to fit the Bergman minimal model with SAAM-II frequently failed to produce acceptable values, and so the study report of effects on insulin sensitivity used Matsuda’s index. We have used the ogttMetrics package to structure the data and compute both Matsuda’s index and the minimal model SI. Figure 4 shows how the diet effects are estimated using these two indices. Confidence intervals are presented for five different contrasts. The left panel of Figure 4 is identical in content to the Insulin sensitivity panel of Figure 3 of Sacks et al. (2014). The right panel shows results based on SI that are qualitatively similar to those found with Matsuda’s index, with the exception of the estimated effect of lowering glycemic index in the context of high overall carbohydrate content. With Matsuda’s index, the 95% confidence interval excludes zero, but this is not observed when SI is used. Further work on optimizing estimation of insulin sensitivity from the 120 minute OGTT protocol is warranted; the empirical Bayes approach of Theodorakis et al. (2017) is of particular interest as individual-level estimation in that procedure borrows strength from information assembled for the cohort as a whole.

Figure 4. Left: Two-sided 95 percent confidence intervals for within-person diet contrasts based on Matsuda’s index.

Right: Analogous confidence intervals for within-person diet contrasts based on SI estimated using ogttMetrics.

Installation and operation

Installation of ogttMetrics can be accomplished using R 3.4 via devtools::install_github("vjcitn/ogttMetrics", dependencies=c("Depends", "Imports", “Suggests")). The key infrastructure components required for ogttMetrics are CRAN package deSolve for minimal model estimation, and Bioconductor package MultiAssayExperiment for data management. The SIexplorer utility employs the shiny package. These key components have extensive dependencies among other CRAN packages, but these dependencies are automatically resolved by the install_github() command given above. All example data analyzed or visualized in this paper are accessible using the data() function. For example, to reproduce Figure 1, use the commands library(ogttMetrics); data(obaSamp); m1 = minmodByID(obaSamp, "1"); plot_OGTT_fit(m1). For Figure 2, in the same session, use QCplots(obasamp).

Users may import data managed in spreadsheets (CSV format) for use with this software. An executable example is available with example(csvImport).