---
title: "Gentamicin (Fuchs 2014)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Gentamicin (Fuchs 2014)}
  %\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::rxode(readModelDb("Fuchs_2014_gentamicin"))$reference`
* Description: `r rxode2::rxode(readModelDb("Fuchs_2014_gentamicin"))$description`
* Article: <https://doi.org/10.1111/bcp.12444>

The packaged model is a two-compartment IV population PK model of
gentamicin in 1449 preterm and term neonates. Body weight enters
allometrically (fixed exponents 0.75 on CL and Q, 1.0 on Vc and Vp),
gestational age (GA) and postnatal age (PNA) enter linearly on CL,
GA enters linearly on Vc, and concomitant dopamine reduces CL by ~12%.

# Population

The model-building cohort comprised 1449 neonates (994 preterm,
455 term) admitted to the Service of Neonatology of Lausanne University
Hospital between December 2006 and October 2011 who received gentamicin
as IV infusion (30 min) for suspected or proven infection, mostly in
combination with amoxicillin. Baseline characteristics (Table 1):
body weight median 2.170 kg (range 0.440-5.510), gestational age median
34 weeks (range 24-42), postnatal age median 1 day (range 0-94),
postmenstrual age median 34.4 weeks; 57.5% male; 20.8% on invasive
ventilation, 59.4% on non-invasive ventilation, 10.6% with patent
ductus arteriosus, 9.4% on dopamine, 1.9% on indomethacin, 0.3% on
furosemide. The conventional initial dose was 3 mg/kg
(December 2006-April 2011) followed by 4 mg/kg (May 2011-October 2011);
plasma drug concentrations were drawn twice (1 h and 12 h after the
first dose) within a routine TDM programme. The model is also useful
for external validation; an independent cohort of 69 neonates was used
in the source paper for accuracy / precision assessment.

The same information is available programmatically via the model's
`population` metadata
(`rxode2::rxode(readModelDb("Fuchs_2014_gentamicin"))$population`).

# Source trace

The per-parameter origin is recorded as an in-file comment next to
each `ini()` entry in
`inst/modeldb/specificDrugs/Fuchs_2014_gentamicin.R`. The table below
collects them in one place for review.

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lcl` (CL)           | 0.089 L/h | Table 2 final model, CL = 0.089 L/h (SE 1%) |
| `lvc` (Vc)           | 0.908 L | Table 2 final model, Vc = 0.908 L (SE 2%) |
| `lq` (Q)             | 0.157 L/h | Table 2 final model, Q = 0.157 L/h (SE 7%) |
| `lvp` (Vp)           | 0.560 L | Table 2 final model, Vp = 0.560 L (SE 4%) |
| `e_wt_cl_q`          | 0.75 (fixed) | Table 2 theta_CL_BW = theta_Q_BW = 0.75 |
| `e_wt_vc_vp`         | 1.0 (fixed) | Table 2 theta_Vc_BW = theta_Vp_BW = 1 |
| `e_ga_cl`            | 1.870 | Table 2 theta_CL_GA = 1.870 (SE 3%) |
| `e_pna_cl`           | 0.054 | Table 2 theta_CL_PNA = 0.054 (SE 6%) |
| `e_conmed_dopa_cl`   | -0.120 | Table 2 theta_CL_DOPA = -0.120 (SE 22%) |
| `e_ga_vc`            | -0.922 | Table 2 theta_Vc_GA = -0.922 (SE 8%) |
| IIV CL (omega^2)     | 0.0755 (CV 28%) | Table 2 BSV CL = 28% CV |
| IIV Vc (omega^2)     | 0.0319 (CV 18%) | Table 2 BSV Vc = 18% CV |
| Cor(CL, Vc)          | 0.87 | Table 2 correlation 87% |
| `addSd`              | 0.10 mg/L | Table 2 additive residual error = 0.10 mg/L (SE 24%) |
| `propSd`             | 0.18 | Table 2 proportional residual error = 18% (SE 1%) |
| Allometric BW form   | `(WT / 2.170)^expt` | Methods "Model-based pharmacokinetic analysis"; reference BW = cohort median 2.170 kg |
| Linear covariate form| `1 + theta * (cov / cov_ref - 1)` | Methods Eq. (2), centred-on-median |
| Categorical covariate form | `1 + theta * IND` | Methods Eq. (5) for categorical covariates |
| Two-compartment IV ODEs | n/a | Results "Base model" + Methods |
| Terminal half-life (reported) | 12.5 h | Discussion paragraph on three-compartment comparison |

# Virtual cohort

Original observed data are not publicly available. The figures below
use virtual subjects whose GA / BW combinations were chosen by the
source paper as illustrative cases (Methods "Model-based
pharmacokinetic analysis", penultimate paragraph): GA = 26, 30, 34, 37,
40 weeks paired with BW = 0.890, 1.080, 2.120, 2.950, 3.580 kg, with
PNA = 1 day (matching the cohort median).

```{r cohort}
set.seed(20140617)

# Five GA / BW pairs from Fuchs 2014 Methods paragraph on simulation
# (Methods "Model-based pharmacokinetic analysis"): the regimens are
# the dosing recommendations from Results paragraph on Figure 4 and
# Table 3 footnote ("Guidelines were as follows..."):
#   GA <= 29 weeks         -> 5   mg/kg q48h
#   30 <= GA <= 34 weeks   -> 4.5 mg/kg q36h
#   GA >= 35 weeks         -> 4   mg/kg q24h
regimen_for <- function(ga) {
  ifelse(ga <= 29, "5 mg/kg q48h",
  ifelse(ga <= 34, "4.5 mg/kg q36h", "4 mg/kg q24h"))
}
dose_per_kg <- function(ga) ifelse(ga <= 29, 5, ifelse(ga <= 34, 4.5, 4))
ii_for     <- function(ga) ifelse(ga <= 29, 48, ifelse(ga <= 34, 36, 24))

representatives <- tibble::tibble(
  ga = c(26, 30, 34, 37, 40),
  bw = c(0.890, 1.080, 2.120, 2.950, 3.580)
) |>
  dplyr::mutate(
    regimen   = regimen_for(ga),
    mg_per_kg = dose_per_kg(ga),
    tau       = ii_for(ga),
    label     = sprintf("GA %2d wk, BW %1.2f kg", ga, bw)
  )
```

```{r cohort-table, echo = FALSE}
representatives |>
  dplyr::select(label, regimen, mg_per_kg, tau) |>
  dplyr::rename(
    "Representative patient" = label,
    "Recommended regimen"    = regimen,
    "Dose (mg/kg)"           = mg_per_kg,
    "Interval (h)"           = tau
  ) |>
  knitr::kable(
  caption = "Representative patients and dosage regimens (Methods + Table 3 footnote)."
)
```

```{r event-table}
# Build event table: each representative gets a 30-min IV infusion
# every tau hours for 3 doses, plus a hourly observation grid spanning
# the full simulation window.
make_subject <- function(idx, row) {
  amt  <- row$bw * row$mg_per_kg
  rate <- amt / 0.5  # 30-minute infusion
  obs_grid <- seq(0, 3 * row$tau, by = 1)
  rxode2::et(amt = amt, rate = rate, ii = row$tau, addl = 2) |>
    rxode2::et(obs_grid) |>
    as.data.frame() |>
    dplyr::mutate(
      id          = idx,
      label       = row$label,
      regimen     = row$regimen,
      WT          = row$bw,
      GA          = row$ga,
      PNA         = 1 / 30.4375,  # median PNA = 1 day expressed in canonical months
      CONMED_DOPA = 0L
    )
}
events <- purrr::map2_dfr(seq_len(nrow(representatives)),
                          split(representatives, seq_len(nrow(representatives))),
                          make_subject)
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
```

# Simulation

```{r simulate}
mod <- rxode2::rxode(readModelDb("Fuchs_2014_gentamicin"))

# Stochastic simulation with population variability (BSV CL 28%, Vc 18%,
# correlation 0.87 and the combined additive/proportional residual error).
# Each representative is replicated 50 times to characterise the prediction
# interval; the resulting 250-subject design renders in well under a minute.
n_rep <- 50L
events_rep <- purrr::map_dfr(seq_len(n_rep), function(rep_i) {
  events |>
    dplyr::mutate(id = id + (rep_i - 1L) * nrow(representatives))
})

sim_rep <- rxode2::rxSolve(mod, events = events_rep,
                           keep = c("label", "regimen", "WT", "GA"))
sim_rep <- as.data.frame(sim_rep)
```

For deterministic typical-value profiles matching Figure 4, zero out
the random effects:

```{r simulate-typical}
mod_typical <- rxode2::zeroRe(mod)
sim_typ <- rxode2::rxSolve(mod_typical, events = events,
                           keep = c("label", "regimen", "WT", "GA"))
sim_typ <- as.data.frame(sim_typ)
```

# Replicate published Figure 4

Figure 4 of Fuchs 2014 shows predicted concentration-time profiles
(typical value plus 90% prediction interval) for the five
representative GA / BW combinations under the dosing regimen chosen
to reach `Cmin <= 1 mg/L` and `Cmax ~ 8 mg/L`.

```{r figure-4}
fig4_band <- sim_rep |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::group_by(label, regimen, time) |>
  dplyr::summarise(
    Q05 = stats::quantile(Cc, 0.05, na.rm = TRUE),
    Q50 = stats::quantile(Cc, 0.50, na.rm = TRUE),
    Q95 = stats::quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(fig4_band, aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.20) +
  geom_line(linewidth = 0.6) +
  geom_hline(yintercept = c(1, 8), linetype = "dashed", colour = "grey50",
             alpha = 0.6) +
  facet_wrap(~ label + regimen, ncol = 2, scales = "free_y") +
  labs(
    x = "Time (h)",
    y = "Gentamicin concentration (mg/L)",
    caption = paste(
      "Replicates Figure 4 of Fuchs 2014.",
      "Solid line: typical-value median.",
      "Shaded ribbon: 5th-95th percentile (90% prediction interval).",
      "Dashed lines: target Cmin <= 1 mg/L and Cmax ~ 8 mg/L."
    )
  )
```

# PKNCA validation

Use PKNCA over the third dosing interval (which approximates steady
state for these neonates given gentamicin's ~12.5 h terminal half-life)
to compute Cmax, Cmin, and the dosing-interval AUC for each
representative regimen. The treatment grouping is the
`representative patient` label so the per-cohort PKNCA results map
back to the published recommendations.

```{r pknca}
# Use the typical-value simulation for the comparison so the results
# correspond to the deterministic prediction shown by the paper's Figure 4
# solid line.
ss_window_for <- function(tau) {
  start_ss <- 2 * tau   # start of third dosing interval
  end_ss   <- 3 * tau
  list(start = start_ss, end = end_ss)
}

# Concentrations: only retain rows within the steady-state window per
# representative subject, and guarantee a row at t = start_ss (pre-dose
# trough) per (id, label).
sim_nca <- sim_typ |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, label, regimen, WT, GA) |>
  dplyr::inner_join(
    representatives |>
      dplyr::mutate(start_ss = 2 * tau, end_ss = 3 * tau) |>
      dplyr::select(label, start_ss, end_ss),
    by = "label"
  ) |>
  dplyr::filter(time >= start_ss & time <= end_ss)

# Defensive: add (id, label, time = start_ss) anchor rows if missing.
anchor_rows <- sim_nca |>
  dplyr::distinct(id, label, regimen, start_ss, end_ss) |>
  dplyr::mutate(time = start_ss, Cc = NA_real_) |>
  dplyr::anti_join(sim_nca |> dplyr::distinct(id, label, time),
                   by = c("id", "label", "time"))
sim_nca <- dplyr::bind_rows(sim_nca, anchor_rows) |>
  dplyr::arrange(label, id, time)

# Need actual Cc value at the anchor; fill from the closest available
# point (typical-value sim is dense so the deviation is negligible).
sim_nca <- sim_nca |>
  dplyr::group_by(id, label) |>
  tidyr::fill(Cc, .direction = "downup") |>
  dplyr::ungroup() |>
  dplyr::filter(!is.na(Cc))

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

# Each representative subject (one id per row of `representatives`) gets a
# single SS-anchor dose at start_ss for PKNCA. The original events table
# uses addl/ii to encode the multi-dose regimen, so the dose rows are
# implicit; we synthesise the SS-anchor dose row directly from the
# representatives metadata.
dose_df <- representatives |>
  dplyr::mutate(
    id    = dplyr::row_number(),
    time  = 2 * tau,
    amt   = bw * mg_per_kg
  ) |>
  dplyr::select(id, time, amt, label)

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

# The intervals data frame carries the grouping column `label` so each
# (start, end) row applies only to the matching subject; without this,
# PKNCA broadcasts every interval to every group and produces duplicate
# / wrong-window rows.
intervals <- representatives |>
  dplyr::transmute(
    label    = label,
    start    = 2 * tau,
    end      = 3 * tau,
    cmax     = TRUE,
    cmin     = TRUE,
    tmax     = TRUE,
    auclast  = TRUE
  ) |>
  as.data.frame()

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

## Terminal half-life check

Fuchs 2014 reports a terminal half-life of 12.5 h (Discussion paragraph
on three-compartment comparison). Computed from the typical-value PK
parameters at the cohort median:

```{r half-life}
typical <- rxode2::rxode(readModelDb("Fuchs_2014_gentamicin"))
ini_vals <- as.data.frame(typical$ini)
get_val <- function(name) ini_vals$est[ini_vals$name == name]
cl_typ <- exp(get_val("lcl"))
vc_typ <- exp(get_val("lvc"))
q_typ  <- exp(get_val("lq"))
vp_typ <- exp(get_val("lvp"))
kel <- cl_typ / vc_typ
k12 <- q_typ  / vc_typ
k21 <- q_typ  / vp_typ
sum_k <- kel + k12 + k21
disc  <- sqrt(sum_k^2 - 4 * kel * k21)
alpha <- (sum_k + disc) / 2
beta  <- (sum_k - disc) / 2
t_half_alpha <- log(2) / alpha
t_half_beta  <- log(2) / beta
cat(sprintf("Typical-value distribution half-life (alpha): %.2f h\n", t_half_alpha))
cat(sprintf("Typical-value terminal half-life (beta):     %.2f h\n", t_half_beta))
cat(sprintf("Fuchs 2014 Discussion reports:                  12.5 h\n"))
```

## Comparison against published Cmax / Cmin / target attainment

Fuchs 2014 chose the per-cohort regimen in Table 3 footnote so the
typical-value Cmax approaches 8 mg/L and Cmin falls below 1 mg/L
(Results "Model validation and simulation"). The table below renders
the simulated Cmax and Cmin for each representative cohort alongside
those targets.

```{r nca_compare}
nca_long <- as.data.frame(nca_res$result) |>
  dplyr::filter(PPTESTCD %in% c("cmax", "cmin")) |>
  dplyr::transmute(
    label_chr = as.character(label),
    param     = PPTESTCD,
    value     = PPORRES
  ) |>
  dplyr::group_by(label_chr, param) |>
  dplyr::summarise(value = mean(value, na.rm = TRUE), .groups = "drop") |>
  tidyr::pivot_wider(names_from = param, values_from = value) |>
  dplyr::rename(label = label_chr)

# Per Figure 4 of Fuchs 2014, the target Cmin <= 1 mg/L and Cmax ~ 8 mg/L
# applies to all five representative regimens (the regimens were chosen
# specifically to hit these targets, Results "Model validation and
# simulation").
published <- representatives |>
  dplyr::transmute(
    label = label,
    regimen = regimen,
    cmax_target = 8.0,
    cmin_target = 1.0
  )

cmp <- nca_long |>
  dplyr::inner_join(published, by = "label") |>
  dplyr::transmute(
    `Representative patient` = label,
    `Regimen`                = regimen,
    `Simulated Cmax (mg/L)`  = round(cmax, 2),
    `Target Cmax (mg/L)`     = cmax_target,
    `Simulated Cmin (mg/L)`  = round(cmin, 2),
    `Target Cmin (mg/L)`     = cmin_target,
    `Cmax within 20% of 8`   = ifelse(abs(cmax - cmax_target) / cmax_target <= 0.20, "yes", "*"),
    `Cmin <= target`         = ifelse(cmin <= cmin_target * 1.20, "yes", "*")
  )

knitr::kable(
  cmp,
  caption = paste(
    "Simulated steady-state Cmax / Cmin at the third dosing interval,",
    "compared with the Cmax ~ 8 mg/L and Cmin <= 1 mg/L targets",
    "underpinning the Table 3 dosing recommendations.",
    "Cells marked '*' deviate by more than 20%."
  ),
  align = c("l", "l", "r", "r", "r", "r", "c", "c")
)
```

Rows marked `*` are expected for the most preterm cohort (GA 26 weeks),
where the q48h interval drives the trough close to (but slightly above)
the 1 mg/L target, and for the larger term cohorts whose typical-value
peak slightly exceeds the 8 mg/L target after three doses. This is
consistent with the source paper's narrative that the chosen regimens
"allow most infants born at a GA above 34 weeks ... to reach target
concentrations" while preterm neonates approach the target boundaries
asymptotically (Results "Model validation and simulation"; Discussion
paragraph on target-attainment narrative). No parameters were tuned to
hit the published targets.

# Assumptions and deviations

- **PNA in canonical months.** Fuchs 2014 reports PNA in days
  (cohort median 1 day, range 0-94). The canonical `PNA` covariate in
  `inst/references/covariate-columns.md` carries months, so the linear
  effect equation is reparameterised inside `model()` as
  `1 + theta * (PNA / pna_ref_months - 1)` with
  `pna_ref_months = 1 / 30.4375` (the cohort median 1 day expressed
  in months). Numerator and denominator carry the same units factor
  and cancel, so the paper's `theta_CL_PNA = 0.054` is unchanged.
- **Allometric exponents wrapped in `fixed()`.** Fuchs 2014 reports
  the allometric body-weight exponents (0.75 on CL / Q, 1.0 on Vc / Vp)
  without standard errors in Table 2 ("-" entries), consistent with
  Methods "Covariate model" which adopts the canonical theoretical
  allometric exponents. They are encoded as `fixed()`.
- **Indomethacin omitted from the final model.** Univariate
  indomethacin co-administration reduced CL by 18%
  (`DeltaOF = -7.8, P = 0.005`), but the bootstrap 95% CI for the
  coefficient included 0 and the paper omitted it from the final model
  (Results "Model validation and simulation"). The library model
  matches the published final structure and does not include an
  indomethacin effect; the indomethacin indicator is documented under
  `covariatesDataExcluded` so its provenance is preserved.
- **Furosemide screened but not retained.** Furosemide
  coadministration reduced CL by 34% in the univariate analysis but
  did not reach the `P < 0.01` threshold used during covariate
  inclusion (`DeltaOF = -6.3, P = 0.012`), and only 5 subjects
  received furosemide (Results "Covariate model" and Discussion). The
  library model does not include a furosemide effect.
- **New `CONMED_DOPA` canonical covariate.** No prior nlmixr2lib
  model registered a dopamine-coadministration indicator, so this
  extraction introduces `CONMED_DOPA` in
  `inst/references/covariate-columns.md` following the established
  `CONMED_<drug>` pattern. Scope is `specific` until a second model
  legitimately reuses it.
- **External validation cohort not simulated.** The 69-subject
  external validation cohort (recruited January-April 2013) is
  described in `population$validation_cohort` but is not used in this
  vignette, which focuses on reproducing the typical-value Cmax / Cmin
  pattern of Figure 4 against the dosing recommendations.
- **No NCA half-life table in the source.** The paper reports a
  single terminal half-life value (12.5 h) qualitatively in the
  Discussion rather than per-cohort NCA half-lives, so the half-life
  comparison above is a single typical-value check rather than a
  per-cohort PKNCA tabulation.
