---
title: "Methods and concepts"
author: "Dieter Menne, Menne Biomed Consulting Tuebingen, dieter.menne@menne-biomed.de"
output: 
  rmarkdown::html_vignette:
    fig_caption: yes
vignette: >
  %\VignetteIndexEntry{Methods and concepts} 
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: breathtestcore.bib
link-citations: true
---


```{r, echo = FALSE, include = FALSE}
suppressPackageStartupMessages(library(knitr))
suppressPackageStartupMessages(library(gridExtra))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(ggplot2))
library(breathtestcore)
knitr::opts_knit$set(unnamed.chunk.label = "btcore_")
opts_chunk$set(comment = NA, fig.width = 4, fig.height = 3)
options(digits = 3)
set.seed(4711)
```


# A heretical view on <sup>13</sup>C breath test for gastric emptying

<sup>13</sup>C breath test data are evaluated in clinical practice and in research as an indirect measure of gastric emptying. The normative literature on the underlying concepts is given in the documentation of `exp_beta()`, `t50_maes_ghoos()` and `t50_bluck_coward()`.  

Here, I present a subjective viewpoint on the pharmakokinetics of breath test time series. After a look at the Maes/Ghoos method (@Ghoos1993, @Maes1998b) as it is currently used,  alternatives are discussed  in the [next section](#secondlook).

## Current concepts

Within the standard framework implemented in this package, data are fitted to an exponential beta function

$${PDR} = m*{dose}*k*\beta(1-e^{-k*minute})^{(\beta - 1)}*e^{-k*minute}$$

Written in R code

```{r, eval = FALSE}
exp_beta = function(minute,dose,m,k,beta) {
     m*dose*k*beta*(1-exp(-k*minute))^(beta-1)*exp(-k*minute)
}
```

The figure below shows examples of emptying curves. The so-called lag time of the Maes/Ghoos method is at the time of the maximum. Calling this "lag" is a highly confusing definition; technically, lag is the time until a first noticeable rise occurs in the curve, best illustrated by the blue curves with `beta = 8`, which has an estimated lag of about 15 minutes. If one looks more closely, even this definition is misleading, because the curves with higher beta are not zero, but behave like a polynomial curve with a power of `beta - 1` at the origin. A real lag would be best illustrated by a curve shifted to the right by `t_lag` and pdr values of 0 for `t < t_lag`, indicating a real physical transport of the meal from the stomach to the small intestine. 

In practice, when <sup>13</sup>C curves are fitted, `beta` is typically in the range of 1.5 and 2.5; values of `beta` near 1 often lead to instabilities and are protected by priors in the Bayesian `breathteststan::stan_fit()` implemented in package `breathteststan`.


```{r, expb, fig.width = 6, echo = FALSE, fig.cap = "Examples of exponential beta functions used to fit breath test curves. _Bottom-up arrow_: Maes/Ghoos lag time. _Top-down arrow_: Maes/Ghoos half emptying time."}
k = c(0.025, 0.015)
beta=c(1.5, 4, 8)
g = expand.grid(minute = 0:300, k = k, beta = beta) %>% 
  mutate(
    pdr = exp_beta(minute, 100, 40, k, beta)
  )
gpar = expand.grid(k = k, beta = beta) %>% 
  mutate(
    t50 = t50_maes_ghoos(.),
    tlag = tlag_maes_ghoos(.),
    y_t50 = exp_beta(t50, 100, 40, k, beta),
    y_tlag = exp_beta(tlag, 100, 40, k, beta)
  )
arrow_1 = arrow(length = unit(0.1, "inches"), type = "closed")
ggplot(g, aes(minute, pdr, color = as.character(beta))) + 
  geom_line() + 
  facet_grid(~paste("k=",k)) + 
  labs(color='beta') + 
  geom_segment(aes(x = t50, y = y_t50 + 5, xend = t50, yend = y_t50 + .3), 
               data = gpar, linewidth = 1, arrow = arrow_1 ) +
  geom_segment(aes(x = tlag, y = y_tlag-10, xend = tlag, yend = y_tlag - .3), 
               data = gpar, linewidth = 1, arrow = arrow_1) 
```

The half-emptying as defined by @Maes1998b is the time where the area under curve (AUC) is half the AUC extrapolated to infinity. In pharmakokinetics terms, it is the time where half of the bioavailable <sup>13</sup>C has been metabolized. This time is determined both by gastric emptying _and_ by the pharmacokinetics of <sup>13</sup>C bound in octanoate/acetate.


```{r, t50, echo = FALSE, fig.cap = "The half-emptying time in the definition of Maes/Ghoos is the time where the area under curve (AUC) is half of the AUC under the curve extrapolated to infinity."}
t50 = (gpar %>% filter(beta == 1.5 & k == 0.015 ))$t50

g %>% filter(beta == 1.5 & k == 0.015 ) %>% 
  ggplot(aes(minute, pdr)) +
  geom_line() +
  geom_area(data = . %>% filter(minute > t50), fill = "lightgray")  +
  geom_area(data = . %>% filter(minute < t50)) 
  
  
```

The AUC _extrapolated to infinity_, which is the denominator in the definition of `t50`, is a vulnerable variable. Even if we accept that the exponential beta function is a reasonable model for breath test time series, extrapolating individual curves from data only slightly longer than `t50` results in ambiguous area estimates, often fails to converge, and even more often gives wildly romantic estimates of `t50`.

To illustrate how brittle the _AUC to infinity_ is, look back at a different function that had been used to fit <sup>13</sup>C time series. @Ghoos1993 used the following Gamma function to fit the data:

$$pdr = a*t^b*e^{-ct}$$
As the figure below shows, this function gives perfectly valid fits of breath test time series. However, the area under curve extrapolated to infinity is infinite, so it cannot be used to compute a half-emptying time from the AUC.

```{r, gamma, fig.cap = "Gamma function from @Ghoos1993 definition as points, and fitted beta exponential as line. Both function look very similar, but the area under the gamma function is infinite and thus cannot be used to compute t50.", echo = FALSE}
gam = function(minute, a, b, c){
  a * minute^b * exp(-c*minute)
}

tibble(minute = seq(0, 300, by = 10)) %>%
  mutate(
    pdr = gam(minute, 1, 1.5, 0.03)
  ) %>% 
  cleanup_data() %>% 
  nls_fit() %>% 
  plot(point_size = 2)


```


## Dieter's soapbox 

I therefore recommend using one of the population method, as implemented in `nlme_fit()` and `stan_fit()` for studies; these methods provide "borrowing strength" to keep the flock of lambs safely together and protect outliers from the big bad wolf.

* If you write or review a paper that uses single curve fits to report study results, rewrite or reject it. Using the function provided in this package is [free as in beer and in speech](https://en.wikipedia.org/wiki/Gratis_versus_libre). The [vignette](https://dmenne.github.io/breathtestcore/articles/data_formats.html) gives details on how to prepare your data with Excel.
* For clinical practice, you have single measurements for which you want an immediate result. You can use the Bayesian method to get regularized estimates; or, if you want to be fancy, prepare a data set of some of your previous cases, add the new patient's data as a mix-in, and use the Mixed-model or the Bayesian fit.

```{r, fig.height = 2.5, fig.cap="The fitting methods available in package `breathtestshiny`", echo = FALSE}
knitr::include_graphics("methods.png")
```


## A second look {#secondlook}

The half-emptying time determined with the <sup>13</sup>C breath test method correlates to some extend with gastric emptying in within-subject measurements, but it is not more than a surrogate for gastric emptying times measured by MRI (with secretion) or scintigraphy (meal-only). For data set `usz_13c_d` in this package, MRI emptying data are available and can be compared with [breath test data](https://dmenne.github.io/breathtestcore/reference/usz_13c_d.html#examples). The type of meal strongly biases the estimates because of the the dependency on lipophilic/lipophobic layering; comparing different meal types is not reliable even within-subject, and less so between-subject.

Many attempts have been made to extract more information from the emptying curve, and to correct for the overall too high values; the Bluck-Coward method implemented in this package is one of them, but the [results](https://dmenne.github.io/breathtestcore/reference/usz_13c_d.html#examples) are are even less consistent with those from MRI than those from the Maes/Ghoos method. @Sanaka2004 mention deconvolution and use an approximate correction for the pharmcological weighting function for the Wagner-Nelson method.

There are many publication that tried an ad-hoc scintigraphic correction (@Keller2009, `t50_maes_ghoos_scintigraphy()`) or one based on pharmacolokinetics (@Bluck2006, `t50_bluck_coward()`). I do not know of any publication which has used rigorous statistical tests such as statistical cross-validation and validation with different meal types to show that some method gives "better" results than the default Maes/Ghoos method.

> In the following, I argue that there is nothing to gain from playing with "better" methods to fit the breath test time series. Independent per-subject information must be available to separate acetate pharmacokinetics, based on concentration gradient diffusion, from gastric empyting, which is physical transport.

### Acetate pharmacology

Assume we do not mix the labelled acetate with the meal, but hide it in acid-resistant coating as it is used for PPI tablets. No <sup>13</sup>C will be recorded as long as  the tablet is in the stomach, but it will quickly release <sup>13</sup>C labelled acetate in the duodenum or small intestine, mimicing a one-time dose. As described by standard pharmacokinetics, this results in a PDR response following a first-order compartmental model. Forgetting about decorative normalization constants, the response is the difference between two exponential functions:

$${foc} =  e^{-k_1 * {minute}} - e^{-k_2 * {minute}}$$
The normalization constant includes clearance and dose and only scales the curve. The correct normalization constants can be found in standard texts on pharmacokinetics (@Gabrielsson2006) or in as a base R function `stats::SSfol()` in @R:PinheiroBates2000. 


```{r, foc, fig.width = 7, fig.asp = 0.4, echo = FALSE, fig.cap="Left: pharmcological PDR response of separated loads released in the smaller bowel after 0, 30 and 90 minutes. Right: sum response when all three loads are released successively."}
foc = function(minute, t_lag, dose, clearance, k1, k2)  {
  norm = dose * k1 * k2/(clearance * (k2-k1)) 
  m_lag = pmax(minute - t_lag, 0)
  f =  norm * ((exp(-k1 * m_lag)) - (exp(-k2 * m_lag)))
}  
ff = expand.grid(minute = 0:300, t_lag_n = c(0, 30, 90)) %>% 
  mutate(
    pdr = foc(minute, t_lag_n, 25, 1, 0.2, 0.05),
    t_lag = as.character(t_lag_n)
  ) 
p1 = ff %>% 
   ggplot(aes(x = minute, y = pdr, group = t_lag, col = t_lag)) +
   geom_line()

p2 = ff %>% 
  group_by(minute) %>% 
  summarize(pdr = sum(pdr)) %>% 
  ggplot(aes(x = minute, y = pdr )) + geom_line()
gridExtra::grid.arrange(p1, p2, ncol = 2, nrow = 1)

```

The left panel in the above plot shows three examples of <sup>13</sup>C fed by an enteric tables, remaining in the stomach for 0, 30 and 90 minutes and then quickly being released in the smaller intestine. Note that we have a real lag caused by physical transport here: before the tablet is in the smaller bowel, there is no response at all. This is different from the polynomial-like pseudo-lag with the beta-exponential function.

The right panel in the above plot shows the summary PDR curve when the three tablets are fed together and by chance are released after 0, 30 and 90 minutes. This mimics the case of a stomach releasing the meal in three chunks of the same size and composition, well known from [pharmacokinetics](https://en.wikipedia.org/wiki/Pharmacokinetics#/media/File:Linear_PK_Example.png) for multiple doses of a drug.

The real flow out of the stomach is not pulse-like, but more like that of a continuous infusion with changing flow. When we assume that the stomach empties with a power-exponential function as in the left panel (see `gastempt::powexp()`), the flow is a wide peak as shown in the right panel below. The PDR response technically is the convolution of the flow with the pulse response shown as red curve (t_lag = 0) in the above plot.

```{r, echo = FALSE}
tempt = 100
beta = 1.8
cap = paste0("Meal volume (left) and flow when gastric emptying function is a power exponential function. Slope was computed with  function `gastempt::powexp_slope()` with inverted sign. 400 ml initial value, tempt = ", tempt, ", beta = ", beta, " Half empyting times t50 are marked by arrows." )
```


```{r, powexp, fig.height = 3, fig.width = 8, echo = FALSE, fig.cap = cap}

# Stolen and simplfied from package gastempt
powexp = function(t, v0, tempt, beta){
  v0 * exp(-(t / tempt) ^ beta)
} 

powexp_slope = function(t, v0, tempt, beta){
  .expr1 <- t/tempt
  .expr4 <- v0 * exp(-.expr1^beta)
  -(.expr4 * (.expr1^(beta - 1) * (beta * (1/tempt))))
}
 
t50_powexp = function(tempt, beta){
  interval = c(0, 1000)
  uniroot(function(t)
        powexp(t, 1, tempt, beta ) - 0.5,
        interval = interval)$root
}

flow = expand.grid(minute = 0:300) %>% 
  mutate(
    flow = - powexp_slope(minute, 400, tempt, beta),
    vol = powexp(minute, 400, tempt, beta)
  ) 


t50_pp = round(t50_powexp(tempt = tempt, beta = beta))
  

p1 = flow %>% 
  ggplot(aes(minute, vol)) + geom_line() +
  ylab("vol in stomach (ml" ) + 
  annotate(geom = "segment", 
    x = t50_pp+50, y = 250, xend = t50_pp, yend = 200,
    linewidth = 1, arrow = arrow_1 )

p2 = flow %>% 
  ggplot(aes(minute, flow)) + geom_line() +
  ylab("Flow out of stomach in ml/m" )

# Convolution
ff_0 = (ff %>%  filter(t_lag == 0))$pdr
cv =  data.frame(minute = 0:300, pdr = 
   convolve(flow$flow, rev(ff_0), type = "open")[1:301]) 

t50_pdr = (cv %>% 
  cleanup_data() %>% 
  nls_fit() %>% 
  coef() %>% 
  filter(method == "maes_ghoos", parameter == "t50"))$value

p3 = cv %>% 
  ggplot(aes(minute, pdr)) + geom_line() +
  ylab("Expected PDF (not normalized)") +
  annotate(geom = "segment",
    x = t50_pdr, y = 45, xend = t50_pdr+50, yend = 40,
    linewidth = 1, arrow = arrow_1 )

gridExtra::grid.arrange(p1, p2, p3, ncol = 3, nrow = 1)


```

The correct half-emptying time `t50` in this example was `r t50_pp` minutes, as determined from the left panel; `t50` determined by fitting with the Maes/Ghoos version in the right panel was `r round(t50_pdr)` minutes. 

The bias from fitting the PDR curve naively is caused by the convolution of mechanical transport in gastric emptying, and pharmacological drug kinetics. The latter is a nuisance parameter and should be eliminated - or is it a parameter of interest in its own?

## How to correct for the pharmacological bias?

> I propose a revised procedure to separate pharmacology from transport. Patient first receive 30 mg acetate in an enteric table together with water. DOB data are recorded for 60 minutes, and after that the normal procedure of meal with 100 mg acetate/octanoate. The resulting PDR will look as follows:

```{r, sim2stage, fig.cap = "Simulated response of a two-stage breath test procedure. The peak within the first hour is from acetate directly released in the small bowel, the second wider peak is the result of the convolution of transport emptying and pharmacokinetics."}
  data.frame(pdr = c(ff_0[0:59]*50, rep(0, 20), cv$pdr)) %>% 
  mutate(
    minute = 0:(n()-1)
  ) %>% 
  filter(minute %% 10 == 0) %>% 
  mutate(
    pdr = pdr + rnorm(n(), 0, 1)
  ) %>% 
  ggplot(aes(x = minute, y = pdr)) +
           geom_point() + geom_line(col = "gray") +
  scale_x_continuous(breaks = seq(0, 360, by = 60))
```

The first narrow peak up to 1 hour is the pharmacological response; by fitting it, the kinetic time constants that broaden the wide peak can be determined, and can be used to remove the bias of the transport component visible in the main peak by mathematical deconvolution.

Can anyone supply a data set for a pilot?

# Why do we need population fits?


## Why no Wagner-Nelson (Sanaka)?

A popular alternative when curve fitting does not work was the Wagner-Nelson (@Sanaka2004) method which uses a non-parametric approach for the initial slope. However, it uses the assumption that the terminal slope is the same for all subjects (k=0.01/min, 0.65/h) which strongly affects the estimate of the half-emptying time, so there is little justification to use the Wagner-Nelson method. 

With some of the functions in this package or with `breathteststan::stan_fit()` in the sister package, all curves can be fitted. I have not seen a single example that fails with the Bayesian method when multiple records are analyzed with a mix-in, so there is no excuse to use the Wagner-Nelson method any longer.

If you must, you can use function [ComputeAndSaveWNFit](https://github.com/dmenne/d13cbreath/blob/master/R/ComputeAndSaveWNFit.R) in my legacy package `D13CBreath`. The function is somewhat awkward to use because it has been written for an application with a tightly integrated database, but feel free to steal the code and run.

## Citations

