---
title: "Gastro-intestinal transit time (GITT) absorption model (Henin 2012)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Gastro-intestinal transit time (GITT) absorption model (Henin 2012)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(nlmixr2lib)
library(PKNCA)
library(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
```

# Model and source

* Citation: Henin E, Bergstrand M, Standing JF, Karlsson MO. A mechanism-based approach for absorption modeling: the gastro-intestinal transit time (GITT) model. The AAPS Journal. 2012;14(2):155-163.
* Article: <https://doi.org/10.1208/s12248-012-9324-y>

Two model files are packaged from this paper because it develops and applies the same GITT framework to two structurally different drugs:

* `modellib("Henin_2012_felodipine")` -- extended-release felodipine with a full drug release and semi-physiological hepatic first-pass structure (Table I).
* `modellib("Henin_2012_diclofenac")` -- enteric-coated diclofenac with a simplified region-specific first-order absorption (Table III).

Both share the tablet transit machinery (paper Equation 1 sigmoid step functions with per-subject inflection points; Table II mean and variance residence times).

# Population

## Felodipine sub-study (Table I)

Cross-over study by Weitschies et al. (2005; Henin 2012 Ref 16) in six healthy volunteers who received magnetically labelled 10 mg extended-release felodipine tablets under fasting and fed conditions. Tablet GI positions and in vivo drug release were monitored by magnetic marker monitoring (MMM), and plasma concentrations were sampled simultaneously. The GITT model was fitted to the felodipine plasma data using MMM residence-time distributions as fixed prior information.

## Diclofenac sub-study (Table III)

Bioequivalence study in 30 healthy adult volunteers who received a single 50 mg oral enteric-coated diclofenac tablet under fasting conditions (unpublished, Rosemont Pharmaceuticals Ltd, UK). Disposition parameters were fixed from an intravenous pediatric study (Korpela 1990; 10 children age 4-6 y receiving 0.5 mg/kg IV) with allometric extrapolation to 70 kg adult weight. The absorption sub-model was estimated on the enteric-coated bioequivalence data. LOQ was 10 ng/mL (= 33.767 nmol/L) and 58% of enteric-coated samples were below LOQ (handled with the M3 method in the paper's original NONMEM fit).

The full population metadata is available programmatically via
`readModelDb("Henin_2012_felodipine")()$population` and
`readModelDb("Henin_2012_diclofenac")()$population`.

# Source trace

The per-parameter origin is recorded as an in-file comment next to each `ini()` entry in the two model files under `inst/modeldb/specificDrugs/`. Highlights:

| Equation / parameter                        | Value                              | Source location                        |
|---------------------------------------------|------------------------------------|----------------------------------------|
| STEP function form                          | STEP(t) = 1/(1+exp(-SIG*(t-IP)))   | Equation 1, page 157 (rewritten from `exp(t*SIG)/(exp(t*SIG)+exp(IP*SIG))` for numerical stability) |
| Sigmoidicity SIG                            | 20 (FIXED; assumption)             | Not reported by Henin 2012; see Assumptions and deviations below |
| Felodipine D1 (fundus)                      | 0.68 mg/h                          | Table I, page 157                      |
| Felodipine D2 = D3 (antrum + PSI)           | 1.91 mg/h                          | Table I, page 157                      |
| Felodipine D4 = D5 (DSI + colon)            | 1.16 mg/h                          | Table I, page 157                      |
| Felodipine K23                              | 0.43 1/h (+5 in SI)                | Table I + footnote a, page 157         |
| Felodipine K34 (fasted / fed)               | 3.48 / 0.81 1/h (+5 in SI)         | Table I + footnote a                   |
| Felodipine K47 = K57                        | 2.87 1/h                           | Table I                                |
| Felodipine K67                              | 1.15 1/h                           | Table I                                |
| Felodipine FA (fasted / fed)                | 0.23 / 0.39                        | Table I                                |
| Felodipine EH                               | 0.50                               | Table I                                |
| Felodipine QH                               | 3.5 L/h/kg^0.75 (FIXED)            | Table I                                |
| Felodipine Vliver                           | 0.0143 L/kg (FIXED)                | Table I                                |
| Felodipine Vcentral / Vperiph1 / Vperiph2   | 20.4 / 88 / 166 L                  | Table I                                |
| Felodipine Qperiph1 / Qperiph2              | 174 / 21.9 L/h                     | Table I                                |
| Felodipine residual error                   | 23% proportional                   | Table I                                |
| Diclofenac CL (per 70kg^0.75)               | 16.5 L/h                           | Table III, page 159                    |
| Diclofenac V1 / Q2 / V2 / Q3 / V3 (per 70kg)| 3.68 / 1.75 / 7.48 / 7.21 / 3.79   | Table III                              |
| Diclofenac KA_PSI / KA_DSI / KA_Col         | 1.06 / 8.64 / 1 (FIXED) 1/h        | Table III                              |
| Diclofenac FA                               | 0.61                               | Table III                              |
| Diclofenac residual error                   | 10 nmol/L additive + 12.6% prop    | Table III                              |
| Fundus MRT (fasted / fed)                   | 0.4 / 1.04 h; VRT 0.46 / 1.09 h^2  | Table II, page 158                     |
| Antrum MRT (fasted / fed)                   | 0.32 / 1.58 h; VRT 0.15 / 2.50 h^2 | Table II                               |
| Proximal SI MRT                             | 1.17 h; VRT 1.37 h^2 (CV 50%)      | Table II                               |
| Distal SI MRT                               | 1.22 h; VRT 1.48 h^2 (CV 58%)      | Table II                               |

# Virtual cohort setup

The per-subject inflection-point covariates (`IP_FA`, `IP_APSI`, `IP_PSI_DSI`, `IP_DSI_C`) are sampled from the paper's Table II log-normal distributions before simulation. The `sample_ips()` helper below packages the sampling and can be reused across both drug scenarios; per SKILL.md discipline it stays inside the vignette (not exported).

```{r helpers}
sample_ips <- function(n, fed = FALSE, seed = 20120728) {
  set.seed(seed)
  # Table II MRT / VRT values (h and h^2)
  mrt_fundus <- if (fed) 1.04 else 0.4
  vrt_fundus <- if (fed) 1.09 else 0.46
  mrt_antrum <- if (fed) 1.58 else 0.32
  vrt_antrum <- if (fed) 2.50 else 0.15
  mrt_psi    <- 1.17;  vrt_psi <- 1.37
  mrt_dsi    <- 1.22;  vrt_dsi <- 1.48

  # log-normal shape from MRT and VRT: eta ~ N(0, VRT); IP = MRT * exp(eta).
  # Table II CV values (100 / 50 / 58) confirm the log-normal parameterisation.
  # For sampling we take sd = sqrt(log(CV^2 + 1)) on the log scale, with CV per
  # region so the median matches MRT.
  cv_to_logsd <- function(cv) sqrt(log(cv^2 + 1))

  sd_fundus <- cv_to_logsd(1.0)   # CV 100%
  sd_antrum <- cv_to_logsd(1.0)   # CV 100%
  sd_psi    <- cv_to_logsd(0.5)   # CV 50%
  sd_dsi    <- cv_to_logsd(0.58)  # CV 58%

  # Residence times per subject (positive, log-normal).
  rt_fundus <- mrt_fundus * exp(stats::rnorm(n, 0, sd_fundus))
  rt_antrum <- mrt_antrum * exp(stats::rnorm(n, 0, sd_antrum))
  rt_psi    <- mrt_psi    * exp(stats::rnorm(n, 0, sd_psi))
  rt_dsi    <- mrt_dsi    * exp(stats::rnorm(n, 0, sd_dsi))

  # Convert per-region residence times to cumulative inflection times.
  ip_fa      <- rt_fundus
  ip_apsi    <- ip_fa   + rt_antrum
  ip_psi_dsi <- ip_apsi + rt_psi
  ip_dsi_c   <- ip_psi_dsi + rt_dsi

  tibble::tibble(
    id         = seq_len(n),
    IP_FA      = ip_fa,
    IP_APSI    = ip_apsi,
    IP_PSI_DSI = ip_psi_dsi,
    IP_DSI_C   = ip_dsi_c
  )
}
```

# Felodipine simulation

The felodipine simulation is a typical-value (`rxode2::zeroRe()`) replication of the extended-release oral dose. To keep the render fast we use a small cohort (`n = 25`) per feeding condition; the primary validation is that the plasma profile is dominated by the region-specific release and absorption events rather than by systemic distribution.

```{r felodipine_events}
mod_fel <- readModelDb("Henin_2012_felodipine")

make_felodipine_cohort <- function(n = 25, fed = 0L, id_offset = 0L, seed = 20120801) {
  ips <- sample_ips(n, fed = as.logical(fed), seed = seed)
  covariates <- ips |>
    dplyr::mutate(id = id + id_offset, WT = 70, FED = fed)

  observe_grid <- seq(0, 24, by = 0.5)

  dose_rows <- covariates |>
    dplyr::transmute(id, time = 0, evid = 1L, amt = 10, cmt = "depot",
                     WT, FED, IP_FA, IP_APSI, IP_PSI_DSI, IP_DSI_C,
                     treatment = ifelse(fed == 1L, "fed", "fasted"))

  obs_rows <- covariates |>
    tidyr::crossing(time = observe_grid) |>
    dplyr::transmute(id, time, evid = 0L, amt = NA_real_, cmt = "central",
                     WT, FED, IP_FA, IP_APSI, IP_PSI_DSI, IP_DSI_C,
                     treatment = ifelse(fed == 1L, "fed", "fasted"))

  dplyr::bind_rows(dose_rows, obs_rows) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

events_fel <- dplyr::bind_rows(
  make_felodipine_cohort(n = 25, fed = 0L, id_offset = 0L,   seed = 20120801),
  make_felodipine_cohort(n = 25, fed = 1L, id_offset = 25L,  seed = 20120802)
)

stopifnot(!anyDuplicated(unique(events_fel[, c("id", "time", "evid")])))
```

```{r felodipine_simulate}
mod_fel_typ <- rxode2::zeroRe(mod_fel)
sim_fel <- rxode2::rxSolve(
  mod_fel_typ, events = events_fel,
  keep = c("treatment", "FED", "WT")
) |> as.data.frame()
```

## Replicates Figure 5 (felodipine VPC by feeding condition)

Figure 5 of the paper shows visual predictive checks of the median felodipine plasma concentration from 1000 simulated replicates for the lag-time, MMM, and GITT approaches. The panel below reproduces the median typical-subject GITT profile for the fasted and fed arms.

```{r felodipine_figure_5}
ggplot(sim_fel, aes(time, Cc, group = id, colour = treatment)) +
  geom_line(alpha = 0.4) +
  facet_wrap(~ treatment) +
  scale_y_continuous(name = "Felodipine plasma concentration (nmol/L)") +
  scale_x_continuous(name = "Time after dose (h)") +
  labs(title = "Figure 5 (Henin 2012) -- typical GITT felodipine profile",
       caption = "Typical-value simulation (zeroRe) with per-subject sampled IPs from Table II.") +
  theme_bw()
```

## PKNCA validation of the felodipine simulation

```{r felodipine_pknca}
sim_nca_fel <- sim_fel |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, treatment)

sim_nca_fel <- dplyr::bind_rows(
  sim_nca_fel,
  sim_nca_fel |> dplyr::distinct(id, treatment) |>
    dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, treatment, time, .keep_all = TRUE) |>
  dplyr::arrange(id, treatment, time)

conc_fel <- PKNCA::PKNCAconc(sim_nca_fel, Cc ~ time | treatment + id)

dose_df_fel <- events_fel |>
  dplyr::filter(evid == 1L) |>
  dplyr::select(id, time, amt, treatment)

dose_fel <- PKNCA::PKNCAdose(dose_df_fel, amt ~ time | treatment + id)

intervals_fel <- data.frame(
  start = 0, end = 24,
  cmax = TRUE, tmax = TRUE, auclast = TRUE, half.life = TRUE
)

nca_fel <- PKNCA::pk.nca(PKNCA::PKNCAdata(conc_fel, dose_fel, intervals = intervals_fel))
```

Summary NCA per feeding condition:

```{r felodipine_nca_summary}
nca_fel_summary <- as.data.frame(nca_fel) |>
  dplyr::group_by(treatment, PPTESTCD) |>
  dplyr::summarise(median = median(PPORRES, na.rm = TRUE), .groups = "drop") |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = median)

nca_fel_summary |>
  dplyr::rename(
    "Feeding" = treatment,
    "Cmax (nmol/L)" = cmax,
    "Tmax (h)" = tmax,
    "AUC0-24 (nmol*h/L)" = auclast,
    "Half-life (h)" = half.life
  ) |>
  knitr::kable(digits = 3, caption = "Median NCA values across simulated felodipine subjects.")
```

Henin 2012 does not report a published NCA table (its validation is graphical VPC), so this table serves as a run-time sanity check rather than a strict comparison. Fed exposure is expected to be higher than fasted because FA is 0.39 (fed) vs 0.23 (fasted).

# Diclofenac simulation

```{r diclofenac_events}
mod_dic <- readModelDb("Henin_2012_diclofenac")

make_diclofenac_cohort <- function(n = 30, id_offset = 0L, seed = 20120901) {
  # Enteric-coated tablet has no drug release in stomach; use only
  # IP_APSI, IP_PSI_DSI, IP_DSI_C. Paper Results reports stomach transit
  # ~ 2 h (range 1.5-3 h) for the diclofenac cohort.
  ips <- sample_ips(n, fed = FALSE, seed = seed) |>
    dplyr::mutate(IP_APSI = pmax(IP_APSI, 1.5),  # enforce reported floor
                  IP_APSI = pmin(IP_APSI, 3.0))  # and ceiling
  covariates <- ips |>
    dplyr::mutate(id = id + id_offset, WT = 70)

  observe_grid <- seq(0, 12, by = 0.25)

  dose_rows <- covariates |>
    dplyr::transmute(id, time = 0, evid = 1L, amt = 50, cmt = "depot",
                     WT, IP_APSI, IP_PSI_DSI, IP_DSI_C,
                     treatment = "enteric-coated 50 mg")

  obs_rows <- covariates |>
    tidyr::crossing(time = observe_grid) |>
    dplyr::transmute(id, time, evid = 0L, amt = NA_real_, cmt = "central",
                     WT, IP_APSI, IP_PSI_DSI, IP_DSI_C,
                     treatment = "enteric-coated 50 mg")

  dplyr::bind_rows(dose_rows, obs_rows) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

events_dic <- make_diclofenac_cohort(n = 30, seed = 20120902)

stopifnot(!anyDuplicated(unique(events_dic[, c("id", "time", "evid")])))
```

```{r diclofenac_simulate}
mod_dic_typ <- rxode2::zeroRe(mod_dic)
sim_dic <- rxode2::rxSolve(
  mod_dic_typ, events = events_dic,
  keep = c("treatment", "WT")
) |> as.data.frame()
```

## Replicates Figure 6 (diclofenac VPC)

Figure 6 of the paper shows the visual predictive checks of diclofenac plasma concentrations for the lag-time, transit-compartment, and GITT models. The panel below reproduces the typical-subject GITT profile.

```{r diclofenac_figure_6}
ggplot(sim_dic, aes(time, Cc, group = id)) +
  geom_line(alpha = 0.4, colour = "#0072B2") +
  scale_y_continuous(name = "Diclofenac plasma concentration (nmol/L)") +
  scale_x_continuous(name = "Time after dose (h)") +
  labs(title = "Figure 6 (Henin 2012) -- typical GITT diclofenac profile",
       caption = "Typical-value simulation (zeroRe) with per-subject sampled IPs from Table II.") +
  theme_bw()
```

## PKNCA validation of the diclofenac simulation

```{r diclofenac_pknca}
sim_nca_dic <- sim_dic |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, treatment)

sim_nca_dic <- dplyr::bind_rows(
  sim_nca_dic,
  sim_nca_dic |> dplyr::distinct(id, treatment) |>
    dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, treatment, time, .keep_all = TRUE) |>
  dplyr::arrange(id, treatment, time)

conc_dic <- PKNCA::PKNCAconc(sim_nca_dic, Cc ~ time | treatment + id)

dose_df_dic <- events_dic |>
  dplyr::filter(evid == 1L) |>
  dplyr::select(id, time, amt, treatment)

dose_dic <- PKNCA::PKNCAdose(dose_df_dic, amt ~ time | treatment + id)

intervals_dic <- data.frame(
  start = 0, end = 12,
  cmax = TRUE, tmax = TRUE, auclast = TRUE, half.life = TRUE
)

nca_dic <- PKNCA::pk.nca(PKNCA::PKNCAdata(conc_dic, dose_dic, intervals = intervals_dic))

nca_dic_summary <- as.data.frame(nca_dic) |>
  dplyr::group_by(treatment, PPTESTCD) |>
  dplyr::summarise(median = median(PPORRES, na.rm = TRUE), .groups = "drop") |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = median)

nca_dic_summary |>
  dplyr::rename(
    "Regimen" = treatment,
    "Cmax (nmol/L)" = cmax,
    "Tmax (h)" = tmax,
    "AUC0-12 (nmol*h/L)" = auclast,
    "Half-life (h)" = half.life
  ) |>
  knitr::kable(digits = 3, caption = "Median NCA values across simulated diclofenac subjects.")
```

Henin 2012 reports the total bioavailability of diclofenac as ~ 61% (Table III FA = 0.61) which is consistent with the literature value cited in the paper (Willis 1979 = ~60%). The simulated AUC scales linearly with FA and the dose.

# Assumptions and deviations

* **STEP-function sigmoidicity SIG is not reported in the paper.** Henin 2012 defines Equation 1 (`STEP(t) = exp(t*SIG)/(exp(t*SIG)+exp(IP*SIG))`, equivalent to the numerically stable `1/(1+exp(-SIG*(t-IP)))` used in the model files) and notes that "the higher the SIG, the steeper the step function", but does not report the numeric value chosen for the fit. `SIG = 20` is used here as a documented assumption; it yields a 10-90 percentile transition width of about 0.22 h around each inflection point. Downstream users can change the value by editing the model file. The paper's supplementary NONMEM control file (which would fix the value unambiguously) is referenced in the paper but is not on disk.
* **Only the 'no return to fundus' subpopulation is encoded.** Henin 2012 fits a 3-component mixture over the antrum-to-fundus return behaviour (subpop 1 = no return, subpop 2 = one return, subpop 3 = two returns) with equal weight (1/3 each) and reports that 7 of 12 felodipine subjects were most likely to belong to the no-return subpopulation. Encoding the mixture would require the nlmixr2 `mixture()` construct which is not exercised elsewhere in this library. This extraction encodes only subpopulation 1 (no return); downstream users who need the fuller mixture behaviour can iterate on the model file to add the two return-to-fundus arms as additional STEP-function pairs (paper Table II shows the fundus and antrum MRTs for each transition).
* **Per-subject IPs are exposed as covariates, not as `ini()` etas.** The paper describes the IPs as sampled from a fixed log-normal distribution `IP = MRT * exp(eta)` with `eta ~ N(0, VRT)`, with the MRT / VRT values themselves fixed to Bergstrand et al. 2009 CPT values (Table II). Because the population parameter distributions are fixed rather than estimated, this extraction exposes each per-subject IP as a covariate (`IP_FA`, `IP_APSI`, `IP_PSI_DSI`, `IP_DSI_C`) and demonstrates the sampling in this vignette via the `sample_ips()` helper. This differs from a typical nlmixr2 formulation that would encode the fixed distributions as `ini()` etas; the covariate form was requested via the initial task sidecar (see task record).
* **Felodipine molecular weight and units.** The paper reports felodipine plasma concentration in the paper as a 23% proportional residual on the observation scale, with LOQ values given implicitly in the tables. This extraction expresses the plasma concentration as `nmol/L` (felodipine molecular weight 384.25 g/mol; conversion to `ng/mL` is `ng/mL = nmol/L * 384.25 / 1000`).
* **Diclofenac absorption units label.** Table III of the paper labels KA_PSI, KA_DSI, KA_Col in `nmol/h`, but dimensionally these are first-order rate constants with units of `1/h`; the `nmol/h` label is a paper typo. The extraction treats them as `h^-1`.
* **Diclofenac disposition scaled from a pediatric IV study.** The disposition parameters in Table III (CL, V1, Q2, V2, Q3, V3) were originally estimated on Ten pediatric IV subjects (Korpela 1990) and allometrically scaled to a 70 kg reference for the adult bioequivalence dataset. The model uses the paper's per-70kg reference values with fixed allometric exponents (0.75 on clearances, 1 on volumes).
* **Diclofenac stomach transit is a single lumped event.** Because the enteric-coated formulation does not release drug in the stomach, IP_FA is not used for the diclofenac model; IP_APSI represents the whole stomach-to-PSI transition (mean ~ 2 h per paper Results). The vignette constrains sampled IP_APSI values to 1.5-3 h to match the range reported in paper Results.
* **Semi-physiological liver blood-plasma ratio.** Table I mentions a fixed blood / plasma concentration ratio of 1.45 for felodipine (referenced from Bergstrand 2009). The extraction operates entirely in plasma-concentration space and does not carry the B / P correction; downstream users who need to translate the hepatic clearance to blood-based clearance can multiply CL_hepatic by 1.45.

# Errata

* The paper mentions an "electronic supplement file" containing an example NONMEM control file; the supplement is not on disk. The extraction is based on the paper text and Tables I, II, III.
* Reference (14) in the paper (Bergstrand et al. 2009, CPT 86:77-83) supplies the Table II MRT / VRT distributions but is likewise not on disk. The Table II values reproduced in the paper are treated as authoritative for the current extraction.
