---
title: "Amoxicillin (Fournier 2018)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Amoxicillin (Fournier 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

```{r meta}
mod_meta <- nlmixr2est::nlmixr(readModelDb("Fournier_2018_amoxicillin"))$meta
```

* Citation: `r mod_meta$reference`
* Description: `r mod_meta$description`
* Article (DOI): <https://doi.org/10.1128/AAC.00505-18>

This vignette validates the packaged `Fournier_2018_amoxicillin` model -- a
two-compartment IV population PK model for amoxicillin in 21 adult ICU burn
patients hospitalized at the CHUV Lausanne Burn Centre -- against the source
publication's Table 1 (baseline demographics), Table 3 (final-model parameter
estimates and covariate equations), Table 4 (PTA stratified by body weight
and renal function), and Figure 2 (visual predictive check).

# Population

The Fournier 2018 analysis prospectively enrolled 21 consecutive adult
patients with severe burns admitted to the Burn Centre of the Centre
Hospitalier Universitaire Vaudois (Lausanne, Switzerland) between October
2013 and October 2016 (ClinicalTrials.gov NCT01965340). Patients received a
course of intravenous amoxicillin alone or in combination with clavulanic
acid for treatment of an infection. The cohort was predominantly male
(76.2%), mean age 50.1 years (SD 24.3, range 16-93). Median body weight on
admission was 72.4 kg (IQR 67.0-83.6; range 60-132). Median total burnt body
surface area was 23% (IQR 12.5-44); 76.2% had inhalation lesions. Mean
SAPS II score was 35.9 (SD 18.6) and median Ryan score was 1 (IQR 1-2).
Median Cockcroft-Gault creatinine clearance was 128 mL/min (IQR 65-150),
with 25% of patients exhibiting augmented renal clearance (CRCL >
150 mL/min). Median ICU length of stay was 23 days (IQR 13.0-39.5) and ICU
mortality was 9.5%. A total of 185 amoxicillin plasma concentrations were
analyzed; 18 of 21 patients contributed a rich kinetic profile (samples at
0, 1, 2, 3, 4, and 5 h after end of infusion) on one occasion.

The same information is available programmatically via the model's
`population` metadata:

```{r population}
str(mod_meta$population)
```

# Source trace

The per-parameter origin is recorded as an in-file comment next to each
`ini()` entry in `inst/modeldb/specificDrugs/Fournier_2018_amoxicillin.R`.
The table below collects them in one place; values come from Fournier 2018
Table 3 "Covariate model mean estimate (RSE, %)" column.

| Parameter / equation                          | Value          | Source location |
|-----------------------------------------------|----------------|-----------------|
| `lcl`  (CL at CRCL = 110 mL/min)              | log(13.6)      | Table 3 row "CL (liter/h)" |
| `lvc`  (V1 at BW = 70 kg)                     | log(9.73)      | Table 3 row "V1 (liters)" |
| `lvp`  (V2)                                   | log(17.6)      | Table 3 row "V2 (liters)" |
| `lq`   (Q)                                    | log(20.1)      | Table 3 row "Q (liters/h)" |
| `e_crcl_cl` (theta_CLCR_CL)                   | 0.57           | Table 3 row "theta_CLCR_CL" |
| `etalcl ~ 0.13028`                            | log(0.373^2+1) | Table 3 row "omega CL (% CV)" = 37.3% |
| `propSd <- 0.37`                              | 0.37           | Table 3 row "Proportional error (%)" = 37 |
| `addSd  <- 0.08`                              | 0.08 mg/L      | Table 3 row "Additive error (mg/liter)" = 0.08 |
| `cl <- exp(lcl + etalcl) * (1 + e_crcl_cl * (CRCL/110 - 1))` | n/a | Table 3 footnote: `TV CL = CL * [1 + theta_CLCR_CL * (CL_CR - 110)]` (interpreted as divisively-centered; see Assumptions) |
| `vc <- exp(lvc) * (WT / 70)`                  | n/a            | Table 3 footnote: `TVV1 = V1 * (BW/70)` |
| `d/dt(central) ... d/dt(peripheral1)`         | n/a            | Results "Population PK model" paragraph 1 ("two-compartment model") |
| `Cc ~ add(addSd) + prop(propSd)`              | n/a            | Results "Population PK model" paragraph 2 ("combined additive and proportional residual error model") |

# Virtual cohort

The original observed amoxicillin concentrations are not publicly available.
The virtual cohort below approximates the published trial demographics and
the Fournier 2018 Methods "Dosing simulations" body-weight distribution:
body weight is log-normally distributed around 89 kg with omega_BW = 0.184
on the log scale (89 * exp(N(0, 0.184^2)); reflects the average BW during
amoxicillin therapy in a related 39-patient burn cohort cited in Methods).
Each subject receives 1 g amoxicillin IV every 8 h infused over 30 min (the
Table 4 reference regimen used for the body-weight/PTA sensitivity
analysis). Six CRCL bands span the renal-function spectrum the paper
simulated (15, 30, 60, 100, 150, 200 mL/min); we use 50 subjects per band
(300 total, well under the 200-per-arm cap) so the band-level VPC remains
comparable across the augmented-CRCL spectrum that motivated the paper's
PTA analysis.

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

# Six CRCL bands matching Fournier 2018 Methods "Dosing simulations".
crcl_bands <- c(15, 30, 60, 100, 150, 200)
n_per_band <- 50L

# Per-band log-normal body weight around 89 kg (paper's Methods
# parameterization: BW = 89 * exp(eta_BW), eta_BW ~ N(0, 0.184^2)),
# truncated to a clinically plausible 40-200 kg range.
make_cohort <- function(crcl, n, id_offset = 0L) {
  wt <- 89 * exp(rnorm(n, mean = 0, sd = 0.184))
  wt <- pmin(pmax(wt, 40), 200)
  tibble::tibble(
    id    = id_offset + seq_len(n),
    band  = sprintf("CRCL = %d", crcl),
    WT_kg = wt,
    CRCL  = crcl
  )
}

cov_tab <- dplyr::bind_rows(lapply(
  seq_along(crcl_bands),
  function(i) make_cohort(
    crcl      = crcl_bands[[i]],
    n         = n_per_band,
    id_offset = (i - 1L) * n_per_band
  )
))

# Dosing: 1 g (= 1000 mg) IV infusion over 30 min, q8h x 30 dose cycle
# (~10 days; matches the Fournier 2018 "steady-state (i.e., after 10 days
# of therapy)" simulation horizon). Observations on dense grid through the
# first dosing interval at steady state for the VPC.
inf_duration_h <- 0.5
dose_amt_mg    <- 1000
dose_interval  <- 8
n_doses        <- 30L
ss_start_h     <- (n_doses - 1L) * dose_interval

dose_times      <- seq(0, by = dose_interval, length.out = n_doses)
obs_grid_within <- seq(0, dose_interval, by = 0.25)
sample_times    <- ss_start_h + obs_grid_within

make_subject <- function(idx, row) {
  doses <- tibble::tibble(
    id   = idx,                 time = dose_times,
    evid = 1L,                  amt  = dose_amt_mg,
    rate = dose_amt_mg / inf_duration_h,
    cmt  = "central",           dv   = NA_real_
  )
  obs <- tibble::tibble(
    id   = idx,                 time = sample_times,
    evid = 0L,                  amt  = NA_real_,
    rate = NA_real_,
    cmt  = "central",           dv   = NA_real_
  )
  dplyr::bind_rows(doses, obs) |>
    dplyr::mutate(
      WT   = row$WT_kg,
      CRCL = row$CRCL,
      band = row$band
    ) |>
    dplyr::arrange(time, dplyr::desc(evid))
}

events <- dplyr::bind_rows(lapply(seq_len(nrow(cov_tab)), function(i) {
  make_subject(idx = cov_tab$id[[i]], row = cov_tab[i, ])
}))

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

# Simulation

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

sim_stoch <- rxode2::rxSolve(
  object = mod, events = events,
  keep   = c("WT", "CRCL", "band")
) |>
  as.data.frame()

sim_typical <- rxode2::rxSolve(
  object = rxode2::zeroRe(mod), events = events,
  keep   = c("WT", "CRCL", "band")
) |>
  as.data.frame()
```

# Replicate published figures

## Figure 2 -- visual predictive check at steady state

The Fournier 2018 Figure 2 prediction-corrected VPC overlays observed
concentrations onto simulated medians and 5th/95th percentiles across the
full study time range (0-24 h within a dosing interval). The simulated
cohort here is grouped by CRCL band so the renal-function effect on
exposure is visible.

```{r figure-2}
# Replicates Figure 2 of Fournier 2018: VPC of amoxicillin concentrations
# at steady state (10 days of 1 g q8h). The published figure shows
# concentrations log-scaled from ~0.1 to >100 mg/L over a dosing interval.
sim_stoch |>
  dplyr::filter(time >= ss_start_h) |>
  dplyr::mutate(t_in_int = time - ss_start_h) |>
  dplyr::group_by(band, t_in_int) |>
  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(aes(t_in_int, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95),
              fill = "gray70", alpha = 0.5) +
  geom_line(linewidth = 0.7) +
  facet_wrap(~band) +
  scale_y_log10() +
  labs(
    x        = "Time within steady-state dosing interval (h)",
    y        = "Amoxicillin in plasma (mg/L, log scale)",
    title    = "Figure 2 -- steady-state VPC by CRCL band",
    subtitle = "Replicates Fournier 2018 Figure 2; 1 g q8h IV over 30 min"
  ) +
  theme_minimal()
```

## Typical-value trajectory at the reference patient

```{r typical-trajectory}
# Reference: CRCL = 100 mL/min, BW = 89 kg cohort (closest to the
# population mean used in the simulations).
ref_band <- "CRCL = 100"
ref_subject_id <- cov_tab |>
  dplyr::filter(band == ref_band) |>
  dplyr::arrange(abs(WT_kg - 70)) |>
  dplyr::slice(1) |>
  dplyr::pull(id)

sim_typical |>
  dplyr::filter(id == ref_subject_id, time >= ss_start_h) |>
  dplyr::mutate(t_in_int = time - ss_start_h) |>
  ggplot(aes(t_in_int, Cc)) +
  geom_line(linewidth = 1) +
  geom_point(size = 1.5) +
  scale_y_log10() +
  labs(
    x        = "Time within steady-state dosing interval (h)",
    y        = "Amoxicillin in plasma (mg/L, log scale)",
    title    = "Typical-value trajectory at the reference patient",
    subtitle = "WT close to 70 kg, CRCL = 100 mL/min; 1 g q8h IV over 30 min"
  ) +
  theme_minimal()
```

# PKNCA validation

PKNCA computes Cmax, Tmax, and AUC over the steady-state dosing interval
on the stochastic cohort, stratified by CRCL band. The Fournier 2018 paper
does not report a tabulated NCA table -- the analysis focuses on PTA at
target f T > MIC -- so the simulated NCA values are reported here as an
internal sanity check that the simulation pipeline produces clinically
plausible exposures across the renal-function spectrum (Cmax in the
~30-150 mg/L range and AUC0-tau decreasing with increasing CRCL).

```{r pknca}
# IMPORTANT: filter on `!is.na(Cc)` only -- adding `time > 0` or `Cc > 0`
# drops the time-zero row PKNCA needs to anchor AUC0-tau.
sim_for_nca <- sim_stoch |>
  dplyr::filter(time >= ss_start_h, !is.na(Cc)) |>
  dplyr::mutate(t_in_int = time - ss_start_h) |>
  dplyr::select(id, t_in_int, Cc, band) |>
  as.data.frame()

# Guarantee a time=0 row per (id, band); pre-dose at steady state ~= Cmin
# but PKNCA needs a non-NA Cc anchor at the start. Use the simulated
# trough Cc at t_in_int = 0 if present (it is, since the obs grid starts
# at 0); otherwise add a synthetic 0 row defensively. The bind_rows /
# distinct pattern keeps an existing 0 row in place.
sim_for_nca <- dplyr::bind_rows(
  sim_for_nca,
  sim_for_nca |>
    dplyr::distinct(id, band) |>
    dplyr::mutate(t_in_int = 0, Cc = 0)
) |>
  dplyr::distinct(id, band, t_in_int, .keep_all = TRUE) |>
  dplyr::arrange(id, band, t_in_int)

conc_obj <- PKNCA::PKNCAconc(
  data    = sim_for_nca,
  formula = Cc ~ t_in_int | band + id,
  concu   = "mg/L",
  timeu   = "hr"
)

# Doses at the start of the steady-state interval (one per subject).
doses_for_nca <- cov_tab |>
  dplyr::mutate(time = 0, amt = dose_amt_mg) |>
  dplyr::select(id, time, amt, band) |>
  as.data.frame()

dose_obj <- PKNCA::PKNCAdose(
  data    = doses_for_nca,
  formula = amt ~ time | band + id,
  doseu   = "mg"
)

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

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

knitr::kable(
  summary(nca_res),
  caption = paste(
    "Simulated steady-state NCA by CRCL band (PKNCA).",
    "1 g q8h IV; AUC0-tau over the 8 h interval."
  )
)
```

## Comparison against published exposures

Fournier 2018 does not tabulate Cmax / AUC by stratum, but the Abstract
and Discussion report typical-patient summaries at CRCL = 110 mL/min and
BW = 70 kg: CL = 13.6 L/h, V1 = 9.7 L, t1/2 = 0.8 h. The typical-value
trajectory above can be verified directly against `mod$meta` and the
simulated `cl`, `vc`, `vp`, `q` from the `sim_typical` output:

```{r model-typical-values}
ref_subject_id_meta <- cov_tab |>
  dplyr::filter(band == "CRCL = 110" | band == "CRCL = 100") |>
  dplyr::arrange(abs(WT_kg - 70)) |>
  dplyr::slice(1) |>
  dplyr::pull(id)

if (length(ref_subject_id_meta) > 0) {
  sim_typical |>
    dplyr::filter(id == ref_subject_id_meta, time >= ss_start_h) |>
    dplyr::summarise(
      cl_Lh    = unique(cl),
      vc_L     = unique(vc),
      vp_L     = unique(vp),
      q_Lh     = unique(q),
      kel_1h   = unique(kel),
      half_life_h = log(2) / unique(kel)
    ) |>
    knitr::kable(caption = paste(
      "Typical PK parameters at the simulated reference subject",
      "(CRCL = 100 mL/min, WT close to 70 kg). Paper Abstract reports",
      "CL = 13.6 L/h, V1 = 9.7 L, t1/2 = 0.8 h at CRCL = 110 mL/min and",
      "BW = 70 kg."
    ))
}
```

As noted in the Assumptions section, the literal Table 3 footnote
equation is replaced by the divisively-centered (median-normalized)
interpretation `TV_CL = 13.6 * (1 + 0.57 * (CRCL/110 - 1))`, which
predicts plausible values across the entire 15-200 mL/min CRCL range
the paper simulated:

* CRCL = 15 mL/min:  `13.6 * (1 + 0.57 * (15/110 - 1))  = 6.91  L/h`
* CRCL = 30 mL/min:  `13.6 * (1 + 0.57 * (30/110 - 1))  = 7.96  L/h`
* CRCL = 60 mL/min:  `13.6 * (1 + 0.57 * (60/110 - 1))  = 10.08 L/h`
* CRCL = 110 mL/min: `13.6 * (1 + 0.57 * 0)             = 13.60 L/h`
* CRCL = 150 mL/min: `13.6 * (1 + 0.57 * (150/110 - 1)) = 16.42 L/h`
* CRCL = 200 mL/min: `13.6 * (1 + 0.57 * (200/110 - 1)) = 19.94 L/h`

These values match the qualitative behaviour described in Fournier 2018
Discussion (clearance rising with augmented renal clearance, falling
with renal impairment) and remain physically positive everywhere; the
literal centered-linear form
`TV_CL = 13.6 * (1 + 0.57 * (CRCL - 110))` would give -374 L/h at
CRCL = 60 and 711 L/h at CRCL = 200, which the paper's PTA simulations
make clear are not what was used.

# Assumptions and deviations

- **CL ~ CRCL covariate equation interpretation -- divisive normalization
  (NOT subtractive centering).** Fournier 2018 Table 3 footnote reports
  `TV_CL = CL * [1 + theta_CLCR_CL * (CL_CR - 110)]` with
  `theta_CLCR_CL = 0.57` as the final-model equation. Taken literally as
  printed (subtractive centering), this yields physically absurd CL
  across the simulated CRCL range:

  * CRCL = 60 mL/min:  `13.6 * (1 + 0.57 * -50) = -374 L/h` (negative)
  * CRCL = 110 mL/min: `13.6 * (1 + 0)          = 13.6 L/h`
  * CRCL = 200 mL/min: `13.6 * (1 + 0.57 * 90)  = 711 L/h`

  Below CRCL ~ 108.2 mL/min the literal form produces negative typical
  CL, contradicting the paper's PTA simulations at CRCL = 15, 30, 60
  mL/min and the Discussion's clinically-sensible CL ratios across the
  augmented-CRCL and renal-impairment strata. The model file uses the
  divisively-centered (median-normalized) interpretation
  `TV_CL = exp(lcl) * (1 + e_crcl_cl * (CRCL/110 - 1))`, which is
  positive everywhere in the simulated 15-200 mL/min CRCL range and
  matches the discussion-level CL pattern (~50% lower at CRCL = 30,
  ~50% higher at CRCL = 200). Both forms agree exactly at CRCL = 110
  (`exp(lcl) = 13.6 L/h`). This follows the precedent set for
  Delattre 2010 amikacin (`inst/modeldb/specificDrugs/Delattre_2010_amikacin.R`)
  where the same absurd-as-printed centered-linear equation was
  operator-resolved to the divisively-centered form. The standing
  policy ("covariate equation absurd as printed -> sensible
  centered/median-normalized interpretation") was applied. No NONMEM
  control stream is available on disk to disambiguate the literal vs.
  median-normalized parameterization, so the resolution depends on the
  only physically tenable reading.

- **V1 ~ BW scaled with allometric exponent 1 (linear), NOT 0.75 / 1.0
  classical allometry.** Fournier 2018 Table 3 footnote reports
  `TVV1 = V1 * (BW/70)` -- a linear scalar with exponent 1. The Methods
  describes this as "allometric scaling" (consistent with the classical
  allometric-scaling convention for volume where the exponent is 1.0 by
  default). No exponent parameter is estimated.

- **CRCL stored under the canonical `CRCL` column despite NOT being
  BSA-normalized.** The canonical `CRCL` column in
  `inst/references/covariate-columns.md` accepts both BSA-normalized
  measured CrCl (mL/min/1.73 m^2) and raw Cockcroft-Gault values
  (mL/min, the source-paper column name `CLCR`). Fournier 2018 uses raw
  Cockcroft-Gault (Methods, "Data collection"). The model documents the
  raw / non-BSA-normalized status in `covariateData[[CRCL]]$units` and
  `notes`. Reference value 110 mL/min (population average) is paper-
  derived; do not compare the magnitude of `e_crcl_cl = 0.57` against
  BSA-normalized reference values listed in the canonical entry -- the
  units differ.

- **No IIV on V1, V2, or Q.** Fournier 2018 Results paragraph 3 reports
  that IIV on V1, V2, and Q "either did not improve the model fit or
  poorly estimated the corresponding random effects" and were not
  retained. The packaged model has IIV only on CL.
  `etalcl ~ log(0.373^2 + 1) = 0.13028` (omega_CL = 37.3% CV converted
  to log-normal variance).

- **Combined additive + proportional residual error.** Methods
  "Population PK model building" describes evaluating proportional,
  additive, and combined error models; Table 3 footnote and Results
  paragraph 3 confirm the combined form was selected. Encoded as
  `Cc ~ add(addSd) + prop(propSd)` with `addSd = 0.08 mg/L` and
  `propSd = 0.37` (Table 3 covariate-model column).

- **Validated CRCL range.** The Fournier 2018 PTA simulations span CRCL
  = 15 mL/min to CRCL = 200 mL/min (Methods, "Dosing simulations"). The
  divisively-centered form used here remains physically positive across
  this entire range. Users are advised to keep simulations within the
  validated 15-200 mL/min range; behaviour at CRCL extremes (< 15 or
  > 250 mL/min) is extrapolation outside the paper's design.

- **Independent (diagonal) OMEGA.** Only one eta is estimated
  (`etalcl`), so the OMEGA matrix is trivially diagonal (1 x 1).

- **Concentration units.** The model uses `mg/L` (paper convention for
  amoxicillin). With dose in `mg` and volumes in `L`, the ratio
  `central / vc` directly gives `mg/L`; no scale factor is applied.

- **Steady-state simulation window.** Fournier 2018 simulated steady
  state "after 10 days of therapy" (Methods, "Dosing simulations"). The
  vignette uses 30 doses of 1 g q8h (10 days) and samples the final
  dosing interval, matching the published simulation horizon.

- **Body-weight distribution.** Methods "Dosing simulations" parameterize
  BW as `89 * exp(eta_BW)` with `eta_BW ~ N(0, 0.184^2)`, representative
  of a 39-patient burn cohort that includes the 21 study patients.
  Encoded directly in the vignette cohort generator.

- **Cohort size.** The vignette uses 50 subjects per CRCL band (300
  total) -- comfortably below the 200-per-arm cap, sufficient for a
  band-stratified VPC and a per-band PKNCA summary. Fournier 2018
  itself simulated 500 patients per condition for its PTA analysis;
  scaling that down is appropriate for a validation vignette whose role
  is to demonstrate the simulation pipeline rather than to re-derive
  PTA.

- **Clavulanic acid not modeled.** The amoxicillin assay used in
  Fournier 2018 did not quantify clavulanic acid; only amoxicillin was
  modeled (Discussion paragraph "Our study had ..."). The packaged model
  reflects this.

- **Race / ethnicity not modeled.** Fournier 2018 does not report race
  composition. The single-centre Swiss tertiary-care setting serves a
  predominantly European population, but the analysis did not test race
  as a covariate and so no race effect is included.
