---
title: "Mecillinam (Koumaki 2023)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Mecillinam (Koumaki 2023)}
  %\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("Koumaki_2023_mecillinam"))$reference`
* Description: `r rxode2::rxode(readModelDb("Koumaki_2023_mecillinam"))$description`
* Article: [Microbiol Spectr 11(2):e03441-22](https://doi.org/10.1128/spectrum.03441-22)
* Underlying concentration data: [Gambertoglio 1980, Antimicrob Agents Chemother 18:952-956](https://doi.org/10.1128/AAC.18.6.952)
* Body-weight reference used by the Monte Carlo simulation: [NCHS Series 3 No. 46 (2021)](https://www.cdc.gov/nchs/data/series/sr_03/sr03-046-508.pdf)

Mecillinam (amdinocillin) is a narrow-spectrum 6-amidinopenicillin with high
activity against *Escherichia coli* and other *Enterobacterales*. Koumaki and
colleagues built a two-compartment population PK model for intravenous
mecillinam and used it in Monte Carlo simulations to propose systemic MIC
breakpoints for intermittent, extended (prolonged), and continuous infusion
regimens.

# Population

The PK parameters were **not** estimated from individual-level data. Koumaki
2023 ("Population pharmacokinetic model") reanalysed the *aggregate*
concentration mean and standard-deviation values published by Gambertoglio
1980 for **12 healthy volunteers** who received a single 10 mg/kg dose of
mecillinam as a 15-min intravenous infusion. The aggregate-data method of
Karakitsios & Dokoumetzidis (2019) was used: at each MCMC step a large virtual
population was drawn from lognormal parameter distributions, the simulated mean
and SD profiles at each time point were computed, and both were fitted to the
observed aggregate mean and SD. Estimation used Stan via RStan 2.19.2 in
R 3.5.1.

Two consequences follow, and both matter when re-using this model:

1. The reported "interindividual variability" absorbs **both** true
   between-subject variability and residual / assay error, because the model
   was fitted to the observed spread of the aggregate data.
2. There is consequently **no separate residual error term** in the published
   model.

The subsequent Monte Carlo simulations were run on 10,000 virtual subjects
(5,000 female, 5,000 male). Because the model is parameterised per kilogram but
doses are given in milligrams, body weight is the dominant source of exposure
variability in those simulations; weights were drawn from lognormal
distributions fitted to the U.S. National Center for Health Statistics
2015-2018 adult weight percentiles (Koumaki 2023 reference 11).

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

# Source trace

The per-parameter origin is recorded as an in-file comment next to each `ini()`
entry in `inst/modeldb/specificDrugs/Koumaki_2023_mecillinam.R`. The table
below collects them in one place for review. Table 1 reports every structural
parameter per kilogram, in mL/min/kg (clearances) and mL/kg (volumes); the
model file stores the same values converted to the model's L/h and L units.

| Equation / parameter | Value as published | Value in `ini()` | Source location |
|---|---|---|---|
| `lcl` | CL = 3.45 mL/min/kg (RSE 1.8%) | `log(0.207)` L/h/kg | Table 1 |
| `lvc` | V1 = 123.12 mL/kg (RSE 16.1%) | `log(0.12312)` L/kg | Table 1 |
| `lq` | Q = 6.74 mL/min/kg (RSE 43.4%) | `log(0.4044)` L/h/kg | Table 1 |
| `lvp` | V2 = 80.56 mL/kg (RSE 21.1%) | `log(0.08056)` L/kg | Table 1 |
| `etalcl` | IIV on CL = 10.2% (RSE 17.1%) | `0.0103503` | Table 1; `log(1 + 0.102^2)` |
| `etalvc` | IIV on V1 = 27.8% (RSE 35.6%) | `0.0744431` | Table 1; `log(1 + 0.278^2)` |
| `etalvp` | IIV on V2 = 36.0% (RSE 35.6%) | `0.1218636` | Table 1; `log(1 + 0.360^2)` |
| IIV on Q | not estimated (blank cell) | *(no eta)* | Table 1; Results: "where no IIV was estimated" |
| `fu` | protein binding assumed 10% | `fixed(0.9)` | Monte Carlo simulations section (reference 12) |
| `propSd` | not estimated | `fixed(0)` | Methods / Table 1: no residual error term |
| Weight scaling | "the model was parameterized per kilogram" | `* WT` (exponent 1) | Monte Carlo simulations section; Table 1 units |
| `d/dt(central)`, `d/dt(peripheral1)` | "a two-compartment PK model" | n/a | "The two-compartment PK model" section |
| `Cu <- fu * Cc` | %*f*T>MIC target on free drug | n/a | Monte Carlo simulations section |

# Virtual cohort

Koumaki 2023 drew body weights from lognormal distributions fitted to the
percentile data of NCHS Series 3 No. 46 (reference 11), stratified by sex. The
paper does not print the resulting mean and SD, so the fit is reproduced here
from the same published percentiles (Tables 3 and 5 of that report, "20 and
over", all race and Hispanic-origin groups). For a lognormal,
`log(x_p) = mu + sigma * qnorm(p)`, so the parameters come from a linear
regression of the log-percentiles on `qnorm(p)`.

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

nchs_p <- c(0.05, 0.10, 0.15, 0.25, 0.50, 0.75, 0.85, 0.90, 0.95)
nchs_female <- c(49.8, 53.9, 57.5, 62.2, 73.1, 88.6, 98.2, 105.3, 119.6)
nchs_male   <- c(61.7, 66.6, 69.9, 75.3, 87.4, 101.9, 110.6, 119.4, 130.3)

fit_lognormal <- function(percentiles, p = nchs_p) {
  z <- stats::qnorm(p)
  co <- stats::coef(stats::lm(log(percentiles) ~ z))
  c(meanlog = unname(co[1]), sdlog = unname(co[2]))
}

ln_female <- fit_lognormal(nchs_female)
ln_male   <- fit_lognormal(nchs_male)

# 200 subjects per arm (100 female + 100 male), the vignette cohort cap.
n_per_sex <- 100L
pop <- dplyr::bind_rows(
  tibble::tibble(SEXF = 1L,
                 WT = stats::rlnorm(n_per_sex, ln_female[["meanlog"]], ln_female[["sdlog"]])),
  tibble::tibble(SEXF = 0L,
                 WT = stats::rlnorm(n_per_sex, ln_male[["meanlog"]], ln_male[["sdlog"]]))
) |>
  dplyr::mutate(id = dplyr::row_number())
```

The fit is checked against the published percentiles it was derived from. The
implied natural-scale means also reproduce the NCHS reported means (77.5 kg for
females, 90.6 kg for males) to better than 0.5%.

```{r cohort-check}
check_fit <- function(par, observed, label) {
  tibble::tibble(
    Sex = label,
    Percentile = paste0(100 * nchs_p, "th"),
    `NCHS (kg)` = observed,
    `Fitted lognormal (kg)` =
      round(stats::qlnorm(nchs_p, par[["meanlog"]], par[["sdlog"]]), 1)
  )
}

dplyr::bind_rows(
  check_fit(ln_female, nchs_female, "Female"),
  check_fit(ln_male,   nchs_male,   "Male")
) |>
  tidyr::pivot_wider(names_from = Percentile,
                     values_from = c(`NCHS (kg)`, `Fitted lognormal (kg)`)) |>
  knitr::kable(caption = "Lognormal body-weight fit vs. the NCHS 2015-2018 percentiles it was derived from.")

tibble::tibble(
  Sex = c("Female", "Male"),
  meanlog = round(c(ln_female[["meanlog"]], ln_male[["meanlog"]]), 4),
  sdlog   = round(c(ln_female[["sdlog"]],   ln_male[["sdlog"]]),   4),
  `Implied mean (kg)` = round(c(
    exp(ln_female[["meanlog"]] + ln_female[["sdlog"]]^2 / 2),
    exp(ln_male[["meanlog"]]   + ln_male[["sdlog"]]^2   / 2)), 1),
  `NCHS reported mean (kg)` = c(77.5, 90.6)
) |>
  knitr::kable(caption = "Fitted lognormal body-weight distributions (NCHS Series 3 No. 46, Tables 3 and 5).")
```

# Simulation

```{r load-model}
mod <- readModelDb("Koumaki_2023_mecillinam")

wt_ref   <- 70
dose_ref <- 10 * wt_ref

# No `id` column: this is a single-subject solve, and supplying `id` would make
# rxode2 treat it as a multi-subject simulation.
ev_typ <- dplyr::bind_rows(
  tibble::tibble(time = 0, amt = dose_ref, dur = 15 / 60, evid = 1L,
                 cmt = "central"),
  tibble::tibble(time = c(seq(0, 0.5, by = 0.005), seq(0.55, 12, by = 0.05)),
                 amt = NA_real_, dur = NA_real_, evid = 0L, cmt = "central")
) |>
  dplyr::mutate(WT = wt_ref, treatment = "10 mg/kg IV, 15-min infusion") |>
  dplyr::arrange(time, dplyr::desc(evid))

# Typical-value (zero-eta) solve. Use the solve-time argument `omega = NA`
# rather than rxode2::zeroRe(): zeroRe() MUTATES state shared with the modeldb
# entry, and whichever kind of solve happens first wins for the rest of the
# session. Calling zeroRe() before the population simulations silently strips
# their between-subject variability; calling it after one leaves the typical
# solve holding a sampled eta set instead of zero. `omega = NA` touches only
# this one solve, so both stay correct regardless of chunk order. The
# assertions below (here and in the PTA chunk) fail the render loudly if
# either property regresses.
sim_typ <- rxode2::rxSolve(mod, events = ev_typ, keep = c("WT", "treatment"),
                           omega = NA) |>
  as.data.frame() |>
  # A single-subject solve returns no `id` column; PKNCA needs it as the
  # subject key, so add it explicitly.
  dplyr::mutate(id = 1L)

stopifnot(
  isTRUE(all.equal(sim_typ$cl[1], 0.207   * wt_ref)),
  isTRUE(all.equal(sim_typ$vc[1], 0.12312 * wt_ref)),
  isTRUE(all.equal(sim_typ$vp[1], 0.08056 * wt_ref))
)
```

## Replicating Figure 1 (single 10 mg/kg dose)

Figure 1 of Koumaki 2023 is a visual-predictive-check-style plot: the mean and
mean +/- SD of the concentrations from simulated profiles, overlaid on the
aggregate mean and SD reported by Gambertoglio 1980 for a single 10 mg/kg
15-min IV infusion. The observed values are shown only as plotted points in the
source figure and are not tabulated anywhere in the paper, so only the model
side can be reproduced here.

A structural note worth making explicit: because every parameter scales
linearly with body weight and the Gambertoglio dose was given per kilogram, the
predicted concentration-time profile for this figure is **independent of body
weight** (`C0 = 10 * WT / (0.12312 * WT)`). Weight only becomes an exposure
driver once fixed milligram doses are given, which is exactly the situation in
the Monte Carlo breakpoint analysis below.

```{r figure-1}
obs_grid <- c(seq(0, 0.5, by = 0.02), seq(0.6, 6, by = 0.1))

ev_single <- pop |>
  dplyr::mutate(amt = 10 * WT, dur = 15 / 60) |>
  tidyr::crossing(tibble::tibble(time = 0)) |>
  dplyr::mutate(evid = 1L, cmt = "central") |>
  dplyr::bind_rows(
    pop |>
      tidyr::crossing(tibble::tibble(time = obs_grid)) |>
      dplyr::mutate(amt = NA_real_, dur = NA_real_, evid = 0L, cmt = "central")
  ) |>
  dplyr::arrange(id, time, dplyr::desc(evid))

sim_single <- rxode2::rxSolve(mod, events = ev_single, keep = c("WT", "SEXF")) |>
  as.data.frame()

sim_single |>
  # Drop the pre-dose t = 0 point (Cc = 0), which has no place on a log axis.
  dplyr::filter(time > 0) |>
  dplyr::group_by(time) |>
  dplyr::summarise(mean_cc = mean(Cc), sd_cc = stats::sd(Cc), .groups = "drop") |>
  ggplot2::ggplot(ggplot2::aes(time, mean_cc)) +
  ggplot2::geom_ribbon(ggplot2::aes(ymin = pmax(mean_cc - sd_cc, 1e-3),
                                    ymax = mean_cc + sd_cc),
                       alpha = 0.2) +
  ggplot2::geom_line(linewidth = 0.9) +
  ggplot2::geom_line(ggplot2::aes(y = pmax(mean_cc - sd_cc, 1e-3)), linetype = "dotted") +
  ggplot2::geom_line(ggplot2::aes(y = mean_cc + sd_cc), linetype = "dotted") +
  ggplot2::scale_y_log10() +
  ggplot2::labs(x = "Time (h)", y = "Mecillinam concentration (mg/L)",
                title = "Figure 1 - mean and mean +/- SD after 10 mg/kg over 15 min",
                caption = "Replicates the model (black / dotted) lines of Figure 1 of Koumaki 2023.")
```

# PKNCA validation

The paper does not print an NCA table, but the Results paragraph gives two
derived quantities that the encoded model must reproduce: a clearance of
**3.45 mL/min/kg** and a steady-state volume of distribution of
**202 mL/kg** ("which can be calculated as the sum of V1 + V2"). Both are
checked here by non-compartmental analysis of a typical-value (no-IIV)
simulation at a 70 kg reference weight, which also verifies the mL/min/kg to
L/h and mL/kg to L unit conversions in `ini()`.

```{r pknca}
# `sim_typ` (the typical-value solve) and `ev_typ` were built in the
# `load-model` chunk, deliberately ahead of any population simulation -- see
# the note there.

# PKNCA input filter: only !is.na(Cc). Do NOT add time > 0 or Cc > 0 -- both
# would drop the time-zero row that anchors AUC0-*.
sim_nca <- sim_typ |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, treatment)

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

dose_df <- ev_typ |>
  dplyr::filter(evid == 1) |>
  dplyr::mutate(id = 1L) |>
  dplyr::select(id, time, amt, dur, treatment)

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | treatment + id,
                             concu = "mg/L", timeu = "h")
# `duration` is required for the IV-infusion parameters: mrt.iv.obs subtracts
# half the infusion duration, and vss.iv.obs = cl.obs * mrt.iv.obs inherits
# that correction. Omitting it silently inflates MRT (and Vss) by dur/2.
dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id, doseu = "mg",
                             duration = "dur")

intervals <- data.frame(
  start = 0, end = Inf,
  cmax = TRUE, tmax = TRUE, aucinf.obs = TRUE, half.life = TRUE,
  cl.obs = TRUE, vss.iv.obs = TRUE
)

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

## Comparison against the published derived values

The reference column converts the paper's per-kilogram values to the 70 kg
reference subject: CL = 3.45 mL/min/kg x 70 kg x 60/1000 = 14.49 L/h, and
Vss = 202 mL/kg x 70 kg / 1000 = 14.14 L.

```{r nca_compare}
published <- tibble::tibble(
  treatment    = "10 mg/kg IV, 15-min infusion",
  cl.obs       = 3.45 * wt_ref * 60 / 1000,
  vss.iv.obs   = 202 * wt_ref / 1000
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated     = nca_res,
  reference     = published,
  by            = "treatment",
  units         = c(cl.obs = "L/h", vss.iv.obs = "L"),
  tolerance_pct = 20
)

knitr::kable(
  cmp,
  caption = "Simulated vs. published derived PK values (70 kg typical subject). * differs from reference by >20%."
)
```

Both agree closely. The residual difference in Vss is a rounding artefact in the
source: V1 + V2 = 123.12 + 80.56 = 203.68 mL/kg, which the Results paragraph
reports as "202 mL/kg". Koumaki 2023 also notes that these values are consistent
with the non-compartmental results of Gambertoglio 1980 (CL 3.5 mL/min/kg,
Vss 230 mL/kg).

```{r derived-extra}
tibble::tibble(
  Quantity = c("CL (mL/min/kg)", "Vss = V1 + V2 (mL/kg)", "Terminal half-life (h)"),
  `Koumaki 2023` = c("3.45", "202 (text) / 203.68 (V1 + V2)", "not reported"),
  `Gambertoglio 1980 (NCA)` = c("3.5", "230", "not reported"),
  `This model (70 kg, typical)` = c(
    sprintf("%.2f", sim_typ$cl[1] / wt_ref * 1000 / 60),
    sprintf("%.2f", (sim_typ$vc[1] + sim_typ$vp[1]) / wt_ref * 1000),
    sprintf("%.3f", as.data.frame(nca_res$result) |>
              dplyr::filter(PPTESTCD == "half.life") |>
              dplyr::pull(PPORRES))
  )
) |>
  knitr::kable(caption = "Derived quantities reported in the text of Koumaki 2023.")
```

The short terminal half-life (about 45 min) is what motivates the paper's
interest in extended and continuous infusion.

# Probability of target attainment

This is the paper's principal result. Target attainment is defined on the
**free** drug (`Cu = fu * Cc`, with 10% protein binding assumed): the fraction
of the dosing interval with free concentration above the MIC must be at least
40% (%*f*T>MIC >= 40%) for intermittent and extended infusion. A regimen is
deemed efficacious when the PTA -- the proportion of subjects attaining the
target -- is at least 90%.

Steady state is evaluated over a dosing interval starting at 48 h. With a
terminal half-life near 45 min this is far beyond steady state (the paper
simulated 5 days), and evaluating an earlier interval changes nothing.

```{r pta-intermittent}
regimens <- tibble::tribble(
  ~regimen,    ~amt, ~ii, ~dur_h,   ~infusion,
  "1 g x 3",   1000,   8, 20 / 60,  "20 min",
  "1 g x 4",   1000,   6, 20 / 60,  "20 min",
  "1.2 g x 4", 1200,   6, 20 / 60,  "20 min",
  "1 g x 3",   1000,   8, 2,        "2 h",
  "1 g x 4",   1000,   6, 2,        "2 h",
  "1.2 g x 4", 1200,   6, 2,        "2 h",
  "1 g x 3",   1000,   8, 4,        "4 h",
  "1 g x 4",   1000,   6, 4,        "4 h",
  "1.2 g x 4", 1200,   6, 4,        "4 h"
)

mics <- c(0.125, 0.25, 0.5, 1, 2, 4, 8, 16)
t_ss <- 48

pta_one_regimen <- function(dose_mg, tau, dur_h) {
  dose_times <- seq(0, t_ss, by = tau)
  obs_times  <- seq(t_ss, t_ss + tau, by = 0.02)
  ev <- dplyr::bind_rows(
    pop |> tidyr::crossing(tibble::tibble(time = dose_times)) |>
      dplyr::mutate(amt = dose_mg, dur = dur_h, evid = 1L, cmt = "central"),
    pop |> tidyr::crossing(tibble::tibble(time = obs_times)) |>
      dplyr::mutate(amt = NA_real_, dur = NA_real_, evid = 0L, cmt = "central")
  ) |>
    dplyr::arrange(id, time, dplyr::desc(evid))

  s <- rxode2::rxSolve(mod, events = ev, keep = c("WT")) |>
    as.data.frame()

  # Guard: the population solve must actually carry between-subject
  # variability. Every parameter scales linearly with WT, so cl / WT is
  # constant across subjects if and only if the etas were dropped (see the
  # note in the load-model chunk). Published IIV on CL is 10.2%.
  per_subject <- unique(data.frame(id = s$id, cl = s$cl, WT = s$WT))
  stopifnot(stats::sd(log(per_subject$cl / per_subject$WT)) > 0.05)

  vapply(mics, function(mic) {
    frac <- tapply(s$Cu > mic, s$id, mean)
    100 * mean(frac >= 0.40)
  }, numeric(1))
}

pta_table2 <- regimens |>
  dplyr::rowwise() |>
  dplyr::mutate(pta = list(pta_one_regimen(amt, ii, dur_h))) |>
  dplyr::ungroup() |>
  tidyr::unnest_wider(pta, names_sep = "_") |>
  dplyr::select(regimen, infusion, dplyr::starts_with("pta_")) |>
  tidyr::pivot_longer(dplyr::starts_with("pta_"), values_to = "PTA") |>
  dplyr::mutate(MIC = rep(mics, times = nrow(regimens))) |>
  dplyr::select(regimen, infusion, MIC, PTA)
```

```{r pta-intermittent-table}
published_t2 <- tibble::tribble(
  ~infusion, ~regimen,    ~MIC,   ~PTA_pub,
  "20 min", "1 g x 3",   0.125,  99.94, "20 min", "1 g x 3",   0.25,  99.66,
  "20 min", "1 g x 3",   0.5,    98.01, "20 min", "1 g x 3",   1,     90.39,
  "20 min", "1 g x 3",   2,      63.99, "20 min", "1 g x 3",   4,     17.93,
  "20 min", "1 g x 3",   8,       0.46, "20 min", "1 g x 3",  16,      0.00,
  "20 min", "1 g x 4",   0.125, 100.00, "20 min", "1 g x 4",   0.25, 100.00,
  "20 min", "1 g x 4",   0.5,    99.98, "20 min", "1 g x 4",   1,     99.75,
  "20 min", "1 g x 4",   2,      96.55, "20 min", "1 g x 4",   4,     72.38,
  "20 min", "1 g x 4",   8,      14.37, "20 min", "1 g x 4",  16,      0.13,
  "20 min", "1.2 g x 4", 0.125, 100.00, "20 min", "1.2 g x 4", 0.25, 100.00,
  "20 min", "1.2 g x 4", 0.5,    99.98, "20 min", "1.2 g x 4", 1,     99.82,
  "20 min", "1.2 g x 4", 2,      97.14, "20 min", "1.2 g x 4", 4,     77.09,
  "20 min", "1.2 g x 4", 8,      22.29, "20 min", "1.2 g x 4",16,      0.57,
  "2 h",    "1 g x 3",   0.125, 100.00, "2 h",    "1 g x 3",   0.25, 100.00,
  "2 h",    "1 g x 3",   0.5,   100.00, "2 h",    "1 g x 3",   1,     99.97,
  "2 h",    "1 g x 3",   2,      98.97, "2 h",    "1 g x 3",   4,     80.72,
  "2 h",    "1 g x 3",   8,      12.90, "2 h",    "1 g x 3",  16,      0.05,
  "2 h",    "1 g x 4",   0.125, 100.00, "2 h",    "1 g x 4",   0.25, 100.00,
  "2 h",    "1 g x 4",   0.5,   100.00, "2 h",    "1 g x 4",   1,    100.00,
  "2 h",    "1 g x 4",   2,     100.00, "2 h",    "1 g x 4",   4,     99.92,
  "2 h",    "1 g x 4",   8,      75.27, "2 h",    "1 g x 4",  16,      4.75,
  "2 h",    "1.2 g x 4", 0.125, 100.00, "2 h",    "1.2 g x 4", 0.25, 100.00,
  "2 h",    "1.2 g x 4", 0.5,   100.00, "2 h",    "1.2 g x 4", 1,    100.00,
  "2 h",    "1.2 g x 4", 2,     100.00, "2 h",    "1.2 g x 4", 4,     99.94,
  "2 h",    "1.2 g x 4", 8,      80.82, "2 h",    "1.2 g x 4",16,     11.49,
  "4 h",    "1 g x 3",   0.125, 100.00, "4 h",    "1 g x 3",   0.25, 100.00,
  "4 h",    "1 g x 3",   0.5,   100.00, "4 h",    "1 g x 3",   1,    100.00,
  "4 h",    "1 g x 3",   2,     100.00, "4 h",    "1 g x 3",   4,     99.88,
  "4 h",    "1 g x 3",   8,      64.98, "4 h",    "1 g x 3",  16,      2.20,
  "4 h",    "1 g x 4",   0.125, 100.00, "4 h",    "1 g x 4",   0.25, 100.00,
  "4 h",    "1 g x 4",   0.5,   100.00, "4 h",    "1 g x 4",   1,    100.00,
  "4 h",    "1 g x 4",   2,     100.00, "4 h",    "1 g x 4",   4,    100.00,
  "4 h",    "1 g x 4",   8,      88.26, "4 h",    "1 g x 4",  16,      9.81,
  "4 h",    "1.2 g x 4", 0.125, 100.00, "4 h",    "1.2 g x 4", 0.25, 100.00,
  "4 h",    "1.2 g x 4", 0.5,   100.00, "4 h",    "1.2 g x 4", 1,    100.00,
  "4 h",    "1.2 g x 4", 2,     100.00, "4 h",    "1.2 g x 4", 4,    100.00,
  "4 h",    "1.2 g x 4", 8,      97.25, "4 h",    "1.2 g x 4",16,     25.61
)

pta_table2 |>
  dplyr::left_join(published_t2, by = c("infusion", "regimen", "MIC")) |>
  dplyr::mutate(cell = sprintf("%.2f / %.1f", PTA_pub, PTA)) |>
  dplyr::select(infusion, MIC, regimen, cell) |>
  tidyr::pivot_wider(names_from = c(infusion, regimen), values_from = cell,
                     names_sep = ", ") |>
  dplyr::rename("MIC (mg/L)" = MIC) |>
  knitr::kable(
    caption = paste(
      "Replicates Table 2 of Koumaki 2023: PTA (%) for intermittent and",
      "prolonged infusion. Each cell is 'published / reproduced'."
    )
  )
```

The reproduction tracks the published table closely across all 72 cells. The
breakpoint conclusions all recover: 1 g x 3 over 20 min supports an MIC
breakpoint of 1 mg/L; 1 g x 4 and 1.2 g x 4 over 20 min support at least
2 mg/L; the 2 h infusions support 2 mg/L (1 g x 3) and 4 mg/L (both QID
regimens); and the 4 h infusions support 4 mg/L (1 g x 3, 1 g x 4) and 8 mg/L
(1.2 g x 4). Cell-level differences of a few percentage points are expected:
this cohort is 200 subjects per arm against the paper's 10,000, so the Monte
Carlo standard error here is roughly 3.5 percentage points near a PTA of 50%.

```{r pta-continuous}
rates_mg_day <- c(150, 500, 1000, 2000, 3500, 4800)

pta_continuous <- function(rate_mg_day) {
  # A single long infusion; the target for continuous regimens is essentially
  # Cu at steady state exceeding the MIC, read at 48 h.
  ev <- dplyr::bind_rows(
    pop |> dplyr::mutate(time = 0, amt = rate_mg_day * 2, dur = 48,
                         evid = 1L, cmt = "central"),
    pop |> dplyr::mutate(time = 48, amt = NA_real_, dur = NA_real_,
                         evid = 0L, cmt = "central")
  ) |>
    dplyr::arrange(id, time, dplyr::desc(evid))

  s <- rxode2::rxSolve(mod, events = ev, keep = c("WT")) |>
    as.data.frame()

  vapply(mics, function(mic) 100 * mean(s$Cu > mic), numeric(1))
}

pta_table3 <- tibble::tibble(rate = rates_mg_day) |>
  dplyr::rowwise() |>
  dplyr::mutate(pta = list(pta_continuous(rate))) |>
  dplyr::ungroup() |>
  tidyr::unnest_wider(pta, names_sep = "_") |>
  tidyr::pivot_longer(dplyr::starts_with("pta_"), values_to = "PTA") |>
  dplyr::mutate(MIC = rep(mics, times = length(rates_mg_day))) |>
  dplyr::select(rate, MIC, PTA)
```

```{r pta-continuous-table}
published_t3 <- tibble::tribble(
  ~rate, ~MIC,  ~PTA_pub,
   150, 0.125, 100.00,  150, 0.25,  98.53,  150, 0.5,  33.93,  150, 1,   0.29,
   150, 2,       0.00,  150, 4,      0.00,  150, 8,     0.00,  150, 16,  0.00,
   500, 0.125, 100.00,  500, 0.25, 100.00,  500, 0.5,  99.91,  500, 1,  64.81,
   500, 2,       2.17,  500, 4,      0.00,  500, 8,     0.00,  500, 16,  0.00,
  1000, 0.125, 100.00, 1000, 0.25, 100.00, 1000, 0.5, 100.00, 1000, 1,  99.91,
  1000, 2,      64.74, 1000, 4,      2.16, 1000, 8,     0.00, 1000, 16,  0.00,
  2000, 0.125, 100.00, 2000, 0.25, 100.00, 2000, 0.5, 100.00, 2000, 1, 100.00,
  2000, 2,      99.91, 2000, 4,     64.81, 2000, 8,     2.17, 2000, 16,  0.00,
  3500, 0.125, 100.00, 3500, 0.25, 100.00, 3500, 0.5, 100.00, 3500, 1, 100.00,
  3500, 2,     100.00, 3500, 4,     99.62, 3500, 8,    50.17, 3500, 16,  0.90,
  4800, 0.125, 100.00, 4800, 0.25, 100.00, 4800, 0.5, 100.00, 4800, 1, 100.00,
  4800, 2,     100.00, 4800, 4,    100.00, 4800, 8,    98.53, 4800, 16, 33.93
)

pta_table3 |>
  dplyr::left_join(published_t3, by = c("rate", "MIC")) |>
  dplyr::mutate(cell = sprintf("%.2f / %.1f", PTA_pub, PTA),
                rate_lab = paste0(format(rate, big.mark = ","), " mg/day")) |>
  dplyr::select(MIC, rate_lab, cell) |>
  tidyr::pivot_wider(names_from = rate_lab, values_from = cell) |>
  dplyr::rename("MIC (mg/L)" = MIC) |>
  knitr::kable(
    caption = paste(
      "Replicates Table 3 of Koumaki 2023: PTA (%) for continuous infusion.",
      "Each cell is 'published / reproduced'."
    )
  )
```

For continuous infusion the free steady-state concentration is
`Cu,ss = fu * Rate / CL`, so the PTA depends only on the **ratio** of the daily
dose to the MIC. The published table shows exactly that symmetry -- 33.93%
appears at both (0.5 mg/L, 150 mg/day) and (16 mg/L, 4,800 mg/day), a 32-fold
change in each -- which is a useful internal check that the underlying model is
linear.

The reproduction matches the published values well for the 500, 1,000, 2,000
and 3,500 mg/day columns, and reproduces the paper's headline conclusions for
2,000 mg/day (MIC 2 mg/L) and 3,500 mg/day (MIC 4 mg/L). It does **not**
reproduce four cells in the two extreme columns, which is examined below.

```{r pta-continuous-diagnosis}
diag <- pta_table3 |>
  dplyr::left_join(published_t3, by = c("rate", "MIC")) |>
  dplyr::mutate(ratio = rate / MIC) |>
  dplyr::filter(PTA_pub > 0.5, PTA_pub < 99.5) |>
  dplyr::group_by(ratio) |>
  dplyr::summarise(Published = mean(PTA_pub), Reproduced = mean(PTA), .groups = "drop") |>
  dplyr::arrange(ratio)

# For continuous infusion PTA = P(CL < fu * Rate / (24 * MIC)) = P(CL < 0.0375 * ratio).
# Anchor a lognormal CL distribution on the two published points at ratio 300
# and 500, then predict the remaining published points.
z  <- stats::qnorm(c(0.3393, 0.6481))
lc <- log(0.0375 * c(300, 500))
sdlog_cl   <- diff(lc) / diff(z)
meanlog_cl <- lc[1] - z[1] * sdlog_cl

diag |>
  dplyr::mutate(
    `CL cut-off (L/h)` = round(0.0375 * ratio, 2),
    `Lognormal fit through the ratio 300 and 500 points` =
      round(100 * stats::plnorm(0.0375 * ratio, meanlog_cl, sdlog_cl), 1)
  ) |>
  dplyr::rename("Dose / MIC ratio" = ratio) |>
  knitr::kable(digits = 2,
               caption = "Continuous-infusion PTA as a function of the dose/MIC ratio.")
```

Because PTA is a monotone function of the dose/MIC ratio, every published cell
must lie on one increasing curve. The published values do not: PTA rises from
2.17% at a ratio of 250 to 33.93% at 300, then only to 50.17% at 437.5 and
64.81% at 500, then jumps to 98.53% at 600. A lognormal clearance distribution
anchored on the ratio-300 and ratio-500 points predicts 24% at ratio 250 (the
paper reports 2.17%) and 75% at ratio 600 (the paper reports 98.53%), so no
smooth clearance distribution can generate the published column. The
reproduction here is smooth and monotone, matches the published values closely
at ratios 437.5 and 500, and agrees that ratios 218.75 and 250 are effectively
unattainable -- these are the ratios contributed by the 500, 1,000, 2,000 and
3,500 mg/day columns. (At the bottom of the range a 200-subject cohort can only
resolve PTA in steps of 0.5%, so the published 0.90% and 2.17% and the
reproduced 0.50% are not meaningfully different.)

The four discrepant cells occur only at ratios 300 and 600, which appear
exclusively in the 150 and 4,800 mg/day columns, and the discrepant values
(33.93, 98.53) are duplicated between those two columns. The most likely
explanation is a transcription error confined to the extreme columns of
Table 3. The practical consequence is that the paper's claim that 4,800 mg/day
covers an MIC of 8 mg/L (published PTA 98.53%) reproduces here at about 85%,
just below the 90% efficacy threshold the paper adopts. **No parameter was
adjusted to reduce this discrepancy.**

# Assumptions and deviations

- **No residual error exists in the published model.** The aggregate-data
  method fitted the observed concentration mean *and* SD, attributing the
  entire observed spread to the lognormal parameter distributions. `propSd` is
  therefore encoded as `fixed(0)` so the model loads and simulates
  deterministically given the etas. The reported IIV consequently absorbs
  both between-subject and residual/assay variability, and anyone re-fitting
  this model on individual-level data must add an estimated residual error term.
- **IIV interpreted as CV%.** Table 1 reports IIV as a percentage of lognormal
  parameter distributions; these are converted to the internal log-scale
  variance with `omega^2 = log(1 + CV^2)`. For the largest value (36% on V2)
  the alternative reading -- treating 0.360 as the log-scale SD directly --
  would give `omega = 0.360` rather than `0.349`, a 3% difference.
- **No IIV on Q.** Table 1 leaves the Q variability cell blank and the Results
  state that "no IIV was estimated" for Q, so no eta was introduced. This is
  the author's structure, not an omission.
- **Per-kilogram scaling with exponent 1.** Table 1 reports every parameter per
  kilogram and the Methods state the model "was parameterized per kilogram", so
  CL, Q, Vc and Vp are multiplied by `WT` with an exponent of 1. No allometric
  (0.75 / 1) exponents were retrofitted.
- **Protein binding is an assumption of the source**, not an estimate:
  "Protein binding was assumed to be 10%", citing the pivmecillinam SmPC. It is
  encoded as `fu <- fixed(0.9)` and drives `Cu`, the free concentration used
  for %*f*T>MIC.
- **Body-weight distribution.** The paper states that lognormal distributions
  were fitted to the NCHS percentile data but does not print the resulting mean
  and SD. The fit was therefore reproduced here from the same published source
  (NCHS Series 3 No. 46, Tables 3 and 5, "20 and over", all race and
  Hispanic-origin groups), which the paper cites as reference 11 but does not
  distribute. The fitted distributions reproduce the published percentiles to
  about 1 kg and the published means to better than 0.5% (see the Virtual
  cohort section), and the close agreement with Table 2 across 72 cells is
  strong evidence that this matches what the authors used. This is a vignette
  simulation input only; no model parameter depends on it.
- **Cohort size.** 200 subjects per arm here versus 10,000 in the paper, per
  this package's vignette cohort cap. The resulting Monte Carlo standard error
  is roughly 3.5 percentage points near a PTA of 50% and smaller in the tails.
- **Steady state evaluated at 48 h** rather than over 5 days. The terminal
  half-life is about 45 minutes, so steady state is reached within the first
  dosing interval and the choice is immaterial.
- **Errata / typographic issues in the source.**
  - Table 3's first MIC row is printed as `1.125`; it is `0.125`. The MIC
    ladder is otherwise identical to Table 2 (0.125 to 16 mg/L), PTA must
    decrease monotonically with MIC, and the dose/MIC scaling symmetry
    confirms it (150 mg/day at MIC 0.125 and 4,800 mg/day at MIC 4 give the
    same ratio of 1,200 and the same PTA of 100%).
  - The Results paragraph reports Vss as "202 mL/kg" but V1 + V2 from Table 1
    is 203.68 mL/kg, a 0.8% rounding discrepancy. The model uses the Table 1
    values.
  - Four cells of Table 3, all in the 150 and 4,800 mg/day columns, are
    inconsistent with the rest of that table and could not be reproduced; see
    the diagnosis in the continuous-infusion section above.
- **Figure 1 observed data are not reproducible.** The Gambertoglio 1980
  aggregate mean and SD values appear only as plotted points in Figure 1 and
  are not tabulated in Koumaki 2023, so only the model curves are replicated.
- No formal erratum or corrigendum to Koumaki 2023 was located.
