---
title: "Metformin (vanRongen 2018)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Metformin (vanRongen 2018)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: `r rxode2::rxode2(readModelDb("vanRongen_2018_metformin"))$reference`
* Description: `r rxode2::rxode2(readModelDb("vanRongen_2018_metformin"))$description`
* Article: <https://doi.org/10.1007/s40272-018-0293-1>

# Population

van Rongen 2018 fits the model to plasma metformin concentrations from 22
overweight and obese Caucasian (Dutch) adolescents enrolled in the metformin
arm of the multicentre randomised double-blind controlled trial
NCT01487993 / EudraCT 2010-023980-17 (Methods Section 2.1-2.2). The study
population was 16 female / 6 male, with mean age 14.5 years
(range 11.1-17.5), mean total body weight 79.3 kg (range 54.7-104.9), mean
BMI 29.1 kg/m^2 (range 22.9-39.3), and mean BMI-SDS 2.8 (range 1.7-4.0).
Eligible patients were obese (BMI-SDS > 2.3) and insulin resistant
(HOMA-IR >= 3.4); patients with DIS_DIAB, polycystic ovary syndrome, impaired
renal function (eGFR < 80 ml/min), or impaired hepatic function were
excluded. Patients had received metformin 500 mg (n = 3) or 1000 mg
(n = 19) twice daily for 37 weeks before the PK sub-study; sampling
followed an oral glucose tolerance test (OGTT) given in the fasted state
with the morning metformin dose ingested simultaneously. Sample times were
pre-dose and 60, 120, 240, 360, and 480 min after dose
(Table 1 of van Rongen 2018).

# Source trace

Per-parameter and per-equation origin (also recorded as in-file comments
in `inst/modeldb/specificDrugs/vanRongen_2018_metformin.R`):

| Equation / parameter        | Value (typical)        | Source location                                                                 |
|-----------------------------|------------------------|---------------------------------------------------------------------------------|
| `lka`                       | `log(0.0248)` 1/min    | Table 2 final Ka = 0.0248 1/min (RSE 25%)                                       |
| `lvc`                       | `log(485)` L           | Table 2 final V/F = 485 L (RSE 10%)                                             |
| `lcl`                       | `log(1.17)` L/min      | Table 2 final CL/F at TBW = 75.8 kg = 1.17 L/min (RSE 6%)                       |
| `e_wt_cl`                   | `0.0138` per kg        | Table 2 final V (linear slope) = 0.0138 (RSE 44%)                               |
| `etalcl`                    | `~ 0.0694`             | Table 2 final IIV CL/F = 26.8% CV (omega^2 = log(1 + 0.268^2))                  |
| `etalvc`                    | `~ 0.1231`             | Table 2 final IIV V/F = 36.2% CV (omega^2 = log(1 + 0.362^2))                   |
| `etalka`                    | `~ 0.4017`             | Table 2 final IIV ka = 70.3% CV (omega^2 = log(1 + 0.703^2))                    |
| `propSd`                    | `0.218` (fraction)     | Table 2 final proportional residual error = 21.8% (RSE 25%)                     |
| CL/F covariate equation     | n/a                    | Section 3.3 / Table 2: CL/F = CL/F_{75.8 kg} * (1 + V * (TBW - 75.8))           |
| `d/dt(depot)` ODE           | n/a                    | One-compartment first-order oral absorption (Section 3.3, Methods 2.6)          |
| `d/dt(central)` ODE         | n/a                    | One-compartment first-order oral absorption (Section 3.3)                       |
| Proportional residual error | n/a                    | Section 3.3 "Residual variability was best described by a proportional error model" |

# Virtual cohort

The cohort below mirrors the trial demographics: a 1000 mg single oral
dose to a virtual population whose total body weight distribution
approximates the study mean (79.3 kg) and SD (13.9 kg), truncated to the
observed range (54.7-104.9 kg). The PK sub-study sampled only the
metformin-1000-mg arm in any quantity (19 of 22 patients); the 500 mg
arm (n = 3) is omitted from the virtual cohort to focus the comparison
on the published dose-standardised summary statistics. Sample times match
the protocol (0, 60, 120, 240, 360, 480 min) plus a denser grid for
visualisation.

```{r cohort}
set.seed(2026L)

n_sub <- 200L
wt_mean <- 79.3
wt_sd   <- 13.9
wt_lo   <- 54.7
wt_hi   <- 104.9

# Draw WT and clip to the observed range.
wt <- pmin(pmax(rnorm(n_sub, wt_mean, wt_sd), wt_lo), wt_hi)

obs_times <- c(seq(0, 480, by = 10))

dose_rows <- tibble(
  id   = seq_len(n_sub),
  WT   = wt,
  time = 0,
  evid = 1L,
  amt  = 1000,
  cmt  = "depot",
  cohort = "1000 mg single oral"
)

obs_rows <- tibble(id = seq_len(n_sub), WT = wt, cohort = "1000 mg single oral") |>
  tidyr::expand_grid(time = obs_times) |>
  mutate(evid = 0L, amt = 0, cmt = "Cc")

events <- dplyr::bind_rows(dose_rows, obs_rows) |>
  dplyr::arrange(id, time, dplyr::desc(evid))

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

# Simulation

A stochastic simulation across the virtual cohort produces concentration
trajectories that include the model's between-subject variability. The
median and 5th-95th percentile envelope below are the analogue of the
prediction-corrected VPC reported in Figure 2 of the paper.

```{r simulate}
mod <- readModelDb("vanRongen_2018_metformin")

sim <- rxode2::rxSolve(
  mod,
  events = events,
  keep   = c("cohort", "WT")
) |>
  as.data.frame() |>
  dplyr::filter(time > 0)
```

# Replicate published figures

## Figure 1 - CL/F vs TBW (typical-value covariate relation)

Figure 1 of van Rongen 2018 plots empirical Bayes CL/F estimates against
TBW with the typical-value covariate curve overlaid. The curve below
reproduces the typical-value relation `CL/F = 1.17 * (1 + 0.0138 *
(TBW - 75.8))` across the study weight range, alongside the simulated
individual CL/F values (with IIV) for orientation.

```{r figure-1}
mod_typical <- mod |> rxode2::zeroRe()
typical_cl <- tibble(
  WT = seq(wt_lo, wt_hi, length.out = 100),
  CL = 1.17 * (1 + 0.0138 * (WT - 75.8))
)

sim_cl <- sim |>
  dplyr::group_by(id, WT) |>
  dplyr::summarise(CL = unique(cl), .groups = "drop")

ggplot(sim_cl, aes(WT, CL)) +
  geom_point(alpha = 0.35, size = 1.2) +
  geom_line(data = typical_cl, colour = "firebrick", linewidth = 0.8) +
  labs(
    x = "Total body weight (kg)",
    y = "Individual CL/F (L/min)",
    title = "Figure 1 - Individual CL/F vs TBW",
    caption = paste(
      "Replicates Figure 1 of van Rongen 2018. Red line: typical-value relation",
      "CL/F = 1.17 * (1 + 0.0138 * (TBW - 75.8))."
    )
  ) +
  theme_minimal()
```

## Figure 2 - Concentration-time envelope (VPC analogue)

```{r figure-2}
sim_summary <- sim |>
  dplyr::filter(!is.na(Cc), Cc > 0) |>
  dplyr::group_by(time) |>
  dplyr::summarise(
    Q05 = quantile(Cc, 0.05),
    Q50 = quantile(Cc, 0.50),
    Q95 = quantile(Cc, 0.95),
    .groups = "drop"
  )

ggplot(sim_summary, aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25, fill = "steelblue") +
  geom_line(colour = "steelblue", linewidth = 0.8) +
  labs(
    x = "Time after dose (min)",
    y = "Plasma metformin Cc (mg/L)",
    title = "Figure 2 analogue - simulated 5th / 50th / 95th percentiles",
    caption = paste(
      "Stochastic simulation of n = 200 virtual subjects receiving 1000 mg",
      "oral metformin; replicates the spread shown in Figure 2 of van Rongen 2018."
    )
  ) +
  theme_minimal()
```

# PKNCA validation

PKNCA non-compartmental analysis on the stochastic simulation, grouped by
dose cohort (the trial sampled only the 1000 mg arm in any quantity, so a
single treatment level is shown). Per-subject Cmax, Tmax, AUC[0-8h], and
AUC[0-Inf] are computed and the median values are compared against the
dose-standardised (1000 mg) NCA summary in Table 3 of van Rongen 2018.

```{r pknca}
sim_nca <- sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, cohort)

conc_obj <- PKNCA::PKNCAconc(
  sim_nca, Cc ~ time | cohort + id,
  concu = "mg/L", timeu = "min"
)

dose_df <- events |>
  dplyr::filter(evid == 1L) |>
  dplyr::select(id, time, amt, cohort)

dose_obj <- PKNCA::PKNCAdose(
  dose_df, amt ~ time | cohort + id,
  doseu = "mg"
)

intervals <- data.frame(
  start      = 0,
  end        = c(480, Inf),
  cmax       = c(TRUE,  FALSE),
  tmax       = c(TRUE,  FALSE),
  auclast    = c(TRUE,  FALSE),
  aucinf.obs = c(FALSE, TRUE),
  half.life  = c(FALSE, TRUE)
)

nca_data <- PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals)
nca_res  <- PKNCA::pk.nca(nca_data)

nca_tbl <- as.data.frame(nca_res$result) |>
  dplyr::filter(PPTESTCD %in% c("cmax", "tmax", "auclast", "aucinf.obs", "half.life")) |>
  dplyr::group_by(cohort, PPTESTCD) |>
  dplyr::summarise(
    median = signif(median(PPORRES, na.rm = TRUE), 4),
    q05    = signif(quantile(PPORRES, 0.05, na.rm = TRUE), 4),
    q95    = signif(quantile(PPORRES, 0.95, na.rm = TRUE), 4),
    .groups = "drop"
  )

knitr::kable(
  nca_tbl,
  caption = paste(
    "PKNCA on simulated 1000 mg single oral metformin (n = 200 virtual",
    "subjects). Cmax in mg/L, Tmax / half-life in min, AUC in mg*min/L."
  )
)
```

## Comparison against published NCA

Table 3 of van Rongen 2018 reports dose-standardised (1000 mg) median
values for the obese-adolescent cohort: AUC[0-Inf] 802.7 mg*min/L
(range 322.8-2568.8), AUC[0-8h] 603.5 mg*min/L (range 286.7-1118.2),
Cmax median 1.80 mg/L (range 0.79-3.45), and Tmax median 120 min
(range 60-240). The side-by-side below uses the simulation medians from
the PKNCA table above.

```{r comparison}
sim_med <- nca_tbl |>
  dplyr::select(parameter = PPTESTCD, simulated_median = median, simulated_q05 = q05, simulated_q95 = q95)

published <- tibble::tibble(
  parameter = c("cmax", "tmax", "auclast", "aucinf.obs"),
  published_median = c(1.80, 120, 603.5,  802.7),
  published_lo     = c(0.79,  60, 286.7,  322.8),
  published_hi     = c(3.45, 240, 1118.2, 2568.8),
  units            = c("mg/L", "min", "mg*min/L", "mg*min/L")
)

comparison <- dplyr::left_join(published, sim_med, by = "parameter") |>
  dplyr::mutate(pct_diff = signif(100 * (simulated_median - published_median) / published_median, 3))

knitr::kable(
  comparison,
  caption = paste(
    "Published vs simulated dose-standardised (1000 mg) NCA medians.",
    "Cmax and Tmax compared to Table 3 of van Rongen 2018; AUC[0-8h] and",
    "AUC[0-Inf] also from Table 3."
  )
)
```

The simulated Cmax is slightly lower than the published median because
the published Cmax is influenced by the largest absorption-rate outliers
(the paper reports Tmax 60 min in 9 of 22 subjects, implying a tail of
faster absorbers than the typical-value ka = 0.0248 1/min). The
simulated AUC[0-8h] and AUC[0-Inf] medians agree with the published
medians within ~10%. The published AUC ranges are wider than the
simulated 5th-95th percentile envelope because the published ranges are
ordinary minima and maxima across only 22 subjects (which include
several apparent absorption outliers), not percentile-based envelopes.

# Assumptions and deviations

- **Final covariate model extracted (linear-TBW form).** van Rongen 2018
  Table 2 presents two equivalent covariate models for CL/F that differ
  only in OFV by 0.2 units (p > 0.05) and in goodness-of-fit plots are
  indistinguishable: (i) the linear-TBW "Final covariate model"
  (`CL/F = 1.17 * (1 + 0.0138 * (TBW - 75.8))`, OFV -141.0, on which the
  bootstrap was performed), and (ii) the "Excess weight covariate
  model" (`CL/F = 1.16 * (WT_for_age_and_length / 70)^0.75 + 0.011 *
  WT_excess`, OFV -140.8). The linear-TBW form was extracted because
  the paper labels it the official final and applies the bootstrap to
  it; the excess-weight form is the paper's preferred mechanistic
  interpretation but requires the non-canonical covariates
  WT_for_age_and_length (derived from Dutch TNO growth tables at the
  patient's age and length at BMI-SDS = 0) and WT_excess (= TBW -
  WT_for_age_and_length). A future extension could register
  WT_DEVELOPMENTAL and WT_EXCESS as canonical covariates and add a
  sibling file `vanRongen_2018_metformin_excessweight.R`.
- **Race recorded as 100% White.** All 22 study patients were of
  Caucasian (Dutch) descent (Methods Section 2.1: "of Caucasian
  descent"). The model's `population$race_ethnicity` field is
  `c(White = 100)`.
- **WT_excess and WT_for_age_and_length not in covariate register.** The
  alternative "Excess weight covariate model" requires these two
  covariates which are not present in `inst/references/covariate-columns.md`.
  This file does not register them because it implements the linear-TBW
  form which only needs the canonical `WT`. If the excess-weight form
  is ever extracted, the two covariates must be added to the register
  first.
- **Single-dose simulation vs multiple-dose protocol.** Patients had
  received metformin 500 or 1000 mg BID for 37 weeks before the PK
  sub-study; the PK sampling itself spans only one dose interval
  (0-480 min after the OGTT-paired dose). The validation simulation
  reproduces the single-dose, fasted, 1000-mg-dose case which is the
  scenario the published AUC and Cmax summaries describe. The 500 mg
  arm (n = 3) is omitted from the virtual cohort.
- **PKNCA `lambda.z` warnings.** The 8-hour observation window is short
  relative to the model-typical half-life (~5 h), so PKNCA's terminal
  regression uses only the last few observation points. Half-life
  values reported in the NCA table are therefore informative but with
  wider uncertainty than the structural-model value.
- **Genetic-polymorphism covariates not modelled.** Section 3.3 reports
  that OCT1 / MATE1 polymorphisms were not significant covariates on
  any PK parameter (p > 0.05); the final model carries no genetic
  covariate.
