---
title: "Tobramycin (Hennig 2008)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Tobramycin (Hennig 2008)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: Hennig S, Norris R, Kirkpatrick CMJ. (2008). Target concentration intervention is needed for tobramycin dosing in paediatric patients with cystic fibrosis -- a population pharmacokinetic study. Br J Clin Pharmacol 65(4):502-510.
* Article: <https://doi.org/10.1111/j.1365-2125.2007.03045.x>

```{r model-meta}
mod_meta <- rxode2::rxode(readModelDb("Hennig_2008_tobramycin"))
mod_meta$description
mod_meta$reference
mod_meta$units
```

# Population

The Hennig 2008 model was developed retrospectively from 318 tobramycin
plasma concentrations measured in 35 paediatric cystic fibrosis (CF)
patients (21 female, 14 male) admitted to Mater Health Services,
Brisbane, Australia between July 2005 and September 2006 (Hennig 2008
Table 1). Mean age was 9.5 years (range 0.5-17.9), mean total body
weight 34.0 kg (6.0-72.6), and mean dose 9.6 mg/kg once-daily IV.
Renal impairment was rare in this cohort -- only one patient had
Cockcroft-Gault creatinine clearance below 50 mL/min, so the model is
most informative for paediatric CF patients with normal serum
creatinine. The retrospective dataset contained on average 9.1 samples
per subject across 4.6 dosing occasions.

The same information is available programmatically via
`readModelDb("Hennig_2008_tobramycin")$population`.

# Source trace

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

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lcl` (log CL/70 kg) | log(6.37) | Hennig 2008 Table 2, Covariate model column |
| `lvc` (log Vc/70 kg) | log(18.70) | Hennig 2008 Table 2, Covariate model column |
| `lq` (log Q/70 kg)   | log(0.39) | Hennig 2008 Table 2, Covariate model column |
| `lvp` (log Vper/70 kg) | log(1.32) | Hennig 2008 Table 2, Covariate model column |
| `e_wt_cl_q` (allometric exponent on CL, Q) | 0.75 | Hennig 2008 Table 2 footnote (f = 3/4 for clearance) |
| `e_wt_vc_vp` (allometric exponent on Vc, Vper) | 1.00 | Hennig 2008 Table 2 footnote (f = 1 for volume) |
| `llag` (log infusion lag) | log(0.40) | Hennig 2008 Table 2, tlag, Covariate model column |
| `ldur` (log infusion duration, fixed) | log(0.5) | Hennig 2008 Table 2, D2 (fixed); Methods page 503 |
| BSV CL (full block element)   | 11.70% CV -> omega^2 = 0.013596 | Hennig 2008 Table 2, Random parameters |
| BSV Vc (full block element)   | 11.66% CV -> omega^2 = 0.013505 | Hennig 2008 Table 2, Random parameters |
| BSV Vper (full block element) | 41.95% CV -> omega^2 = 0.162130 | Hennig 2008 Table 2, Random parameters |
| Cor(CL, Vc)   | 0.73 | Hennig 2008 Table 2, Correlations between BSVs |
| Cor(CL, Vper) | 0.49 | Hennig 2008 Table 2, Correlations between BSVs |
| Cor(Vc, Vper) | 0.27 | Hennig 2008 Table 2, Correlations between BSVs |
| `propSd` (proportional residual SD) | 0.190 | Hennig 2008 Table 2, Residual variability (19.0% CV exponential -> proportional) |
| Two-compartment IV ODE structure | n/a | Hennig 2008 Methods, Population analysis section, page 503 |
| Allometric scaling formula | `(WT/70)^f` with f = 3/4 (CL, Q) and f = 1 (Vc, Vper) | Hennig 2008 Table 2 footnote |

# Virtual cohort

The original observed concentrations are not publicly available. The
figures below use a virtual paediatric CF cohort whose body-weight
distribution approximates Table 1 of Hennig 2008 (mean 34.0 kg, range
6.0-72.6 kg). Once-daily IV infusion doses are computed from a per
mg/kg regimen (the dosing strategies compared in Hennig 2008 Table 3).

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

n_subjects <- 200L

# Approximate the Hennig 2008 Table 1 weight distribution: mean 34 kg,
# range 6-72.6 kg. Use a truncated normal-ish draw to span the observed
# range; the simulation only consumes WT for allometric scaling.
sample_wt <- function(n) {
  wt <- rnorm(n, mean = 34, sd = 16)
  wt <- pmin(pmax(wt, 6.0), 72.6)
  wt
}

# Dosing regimens evaluated in Table 3 of Hennig 2008. Per-kg doses are
# computed from the simulated WT; the absolute 350 mg ("one dose fits all")
# regimen is applied without weight scaling.
regimens <- tibble::tribble(
  ~regimen,        ~dose_mg_per_kg, ~fixed_dose_mg,
  "350 mg",        NA_real_,        350,
  "7.5 mg/kg",     7.5,             NA_real_,
  "10 mg/kg",      10,              NA_real_
)

make_cohort <- function(n, regimen, dose_mg_per_kg, fixed_dose_mg, id_offset = 0L) {
  wt <- sample_wt(n)
  amt <- if (is.na(fixed_dose_mg)) wt * dose_mg_per_kg else rep(fixed_dose_mg, n)
  ids <- id_offset + seq_len(n)

  dose_rows <- tibble::tibble(
    id   = ids,
    time = 0,
    evid = 1L,
    cmt  = "central",
    amt  = amt,
    WT   = wt,
    regimen = regimen
  )

  obs_grid <- seq(0, 24, by = 0.25)
  obs_rows <- tidyr::expand_grid(id = ids, time = obs_grid) |>
    dplyr::left_join(dplyr::select(dose_rows, id, WT, regimen), by = "id") |>
    dplyr::mutate(evid = 0L, cmt = "central", amt = 0)

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

events <- dplyr::bind_rows(
  make_cohort(n_subjects, "350 mg",    NA_real_, 350,      id_offset =   0L),
  make_cohort(n_subjects, "7.5 mg/kg", 7.5,      NA_real_, id_offset = 1000L),
  make_cohort(n_subjects, "10 mg/kg",  10,       NA_real_, id_offset = 2000L)
)

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

n_doses <- sum(events$evid == 1L)
n_obs   <- sum(events$evid == 0L)
cat("Dosing rows:", n_doses, " | Observation rows:", n_obs, "\n")
```

# Simulation

```{r simulate}
mod <- readModelDb("Hennig_2008_tobramycin")
sim <- rxode2::rxSolve(
  mod,
  events = events,
  keep   = c("regimen", "WT")
) |>
  as.data.frame()
```

# Replicate published figures

## Figure 3 -- visual predictive check by dose regimen

Hennig 2008 Figure 3 shows the model's VPC at the 350 mg regimen
(observed grey points and grey 50th-percentile line vs. simulated
median and 5th/95th-percentile envelopes). The chunk below renders the
analogous percentile envelope for each of the three regimens.

```{r figure-3}
vpc_summary <- sim |>
  dplyr::group_by(regimen, time) |>
  dplyr::summarise(
    Q05 = quantile(Cc, 0.05, na.rm = TRUE),
    Q50 = quantile(Cc, 0.50, na.rm = TRUE),
    Q95 = quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(vpc_summary, aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95, fill = regimen), alpha = 0.2) +
  geom_line(aes(colour = regimen)) +
  facet_wrap(~regimen) +
  scale_y_log10() +
  labs(
    x = "Time after dose (h)",
    y = "Tobramycin concentration (mg/L)",
    title = "Replicates Figure 3 of Hennig 2008",
    caption = "Simulated 5th/50th/95th-percentile envelopes by dosing regimen."
  ) +
  theme(legend.position = "none")
```

# PKNCA validation

Compute Cmax and AUC0-24 (the dosing interval) per subject with
PKNCA. Treatment grouping is by `regimen`.

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

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

conc_obj <- PKNCA::PKNCAconc(
  sim_nca, Cc ~ time | regimen + id,
  concu = "mg/L", timeu = "h"
)
dose_obj <- PKNCA::PKNCAdose(
  dose_df, amt ~ time | regimen + id,
  doseu = "mg"
)

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

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

## Comparison against Hennig 2008 Table 3 simulation summaries

Hennig 2008 Table 3 reports, for each dosing regimen, the percentage
of patients whose AUC0-24 falls within the 80-125 mg*h/L target range
and whose Cmax falls within the 24-38 mg/L target range, based on
1000 Monte Carlo simulations. The chunk below compares the simulated
percentages from this packaged model with the published values.

```{r table-3-comparison}
res_tbl <- as.data.frame(nca_res$result)

per_subject <- res_tbl |>
  dplyr::filter(PPTESTCD %in% c("auclast", "cmax")) |>
  tidyr::pivot_wider(
    id_cols     = c("regimen", "id"),
    names_from  = PPTESTCD,
    values_from = PPORRES
  )

simulated_pct <- per_subject |>
  dplyr::group_by(regimen) |>
  dplyr::summarise(
    pct_auc_in_target  = mean(auclast >= 80   & auclast <= 125, na.rm = TRUE) * 100,
    pct_auc_below_80   = mean(auclast < 80,                     na.rm = TRUE) * 100,
    pct_auc_above_125  = mean(auclast > 125,                    na.rm = TRUE) * 100,
    pct_cmax_in_target = mean(cmax    >= 24   & cmax    <= 38,  na.rm = TRUE) * 100,
    .groups = "drop"
  )

published_pct <- tibble::tribble(
  ~regimen,    ~pct_auc_in_target_pub, ~pct_auc_below_80_pub, ~pct_auc_above_125_pub, ~pct_cmax_in_target_pub,
  "350 mg",    37.9,                    36.6,                  25.5,                   29.8,
  "7.5 mg/kg", 21.7,                    78.1,                   0.2,                   17.3,
  "10 mg/kg",  72.1,                    22.2,                   5.7,                   91.6
)

comparison <- published_pct |>
  dplyr::left_join(simulated_pct, by = "regimen") |>
  dplyr::select(regimen,
                pct_auc_in_target_pub,  pct_auc_in_target,
                pct_auc_below_80_pub,   pct_auc_below_80,
                pct_auc_above_125_pub,  pct_auc_above_125,
                pct_cmax_in_target_pub, pct_cmax_in_target)

knitr::kable(
  comparison,
  digits  = 1,
  caption = paste("Comparison of simulated AUC0-24 and Cmax target-attainment percentages",
                  "with Hennig 2008 Table 3. The 'pub' columns are the published values; the",
                  "remaining columns are computed from the packaged model with",
                  paste0("n = ", n_subjects, " subjects per regimen."))
)
```

The simulated percentages should track the published Table 3 values
within reasonable Monte Carlo noise. AUC target attainment matches the
published values across all three regimens (within ~10 percentage
points). Cmax target attainment matches well for the 350 mg flat dose
but diverges for the per-kg regimens; in particular, the 7.5 mg/kg
arm's published 17.3% in the 24-38 mg/L target band differs from the
simulated value by more than 20 percentage points. The most likely
explanation is a difference in how the paper's simulated "Cmax" was
defined (e.g., peak observed at a specific end-of-infusion sampling
time rather than the maximum over the dosing interval). Because Cmax
under per-kg dosing in a fully-allometric model is approximately
weight-independent, with the typical value
(7.5 mg/kg x 70 kg / Vc_70 = 28 mg/L) sitting near the centre of the
24-38 mg/L band, almost any draw from the BSV distribution lands in
target -- which is what the simulation shows. The model parameters and
structure are reproduced faithfully from the source; per the
extraction skill's policy, parameters are not tuned to chase a
validation metric. Users who need to match the paper's Cmax-distribution
results closely should re-derive Cmax according to the paper's exact
sampling convention before comparing.

# Assumptions and deviations

- **Between-occasion variability (BOV) on CL excluded.** Hennig 2008
  estimated BOV on CL of 6.47% CV in addition to BSV (11.70% CV).
  Following the precedent set by `Faelens_2021_infliximab` -- which
  also drops a low-magnitude per-occasion variance for nlmixr2lib
  portability -- the packaged model omits BOV. The omitted variance is
  small relative to BSV, so simulated population-level statistics
  closely match the published values. Re-introducing BOV would require
  a per-occasion eta multiplexed by an occasion indicator (cf.
  `Wilkins_2008_rifampicin` for that pattern).
- **Body-weight distribution.** The vignette draws WT from a truncated
  normal with mean 34 kg matching the published baseline mean (Table 1)
  and clamped to the 6-72.6 kg study range. The original individual WT
  values are not published.
- **Renal function.** Only one patient in the source cohort had reduced
  creatinine clearance, and serum creatinine / CrCL were tested as
  covariates but not retained. The packaged model therefore does not
  include any renal covariate; users with renally impaired patients
  (e.g. eGFR substantially below the paediatric CF norm) should
  interpret simulations with care.
- **Allometric exponents fixed.** The Hennig 2008 Table 2 footnote
  states "f = 3/4 for clearance, f = 1 for volume" without confidence
  intervals, indicating these were fixed. The packaged model wraps
  `e_wt_cl_q` and `e_wt_vc_vp` in `fixed()` to mirror this.
- **Infusion duration fixed at 30 min.** Hennig 2008 fixed D2 = 0.5 h
  per the hospital protocol; the packaged model imposes
  `dur(central) <- 0.5` h so users do not need to supply a `rate` or
  `dur` column on each dose record.
