Skip to contents

Mecillinam (Koumaki 2023)

Source: vignettes/articles/Koumaki_2023_mecillinam.Rmd
Koumaki_2023_mecillinam.Rmd

Model and source

  • Citation: Koumaki V, Dokoumetzidis A, Angelerou MGF, Baka S, Balakrishnan I, Tsakris A. (2023). Pharmacokinetic/pharmacodynamic determination of systemic MIC breakpoints for intermittent, extended, and continuous infusion dosage regimens of mecillinam. Microbiology Spectrum 11(2):e03441-22. doi:10.1128/spectrum.03441-22. The underlying concentration-time data (aggregate mean and SD, 12 healthy volunteers, single 10 mg/kg 15-min IV infusion) are from Gambertoglio JG, Barriere SL, Lin ET, Conte JE Jr. (1980). Pharmacokinetics of mecillinam in healthy subjects. Antimicrob Agents Chemother 18:952-956. doi:10.1128/AAC.18.6.952 (Koumaki 2023 reference 8). The aggregate-data estimation methodology is Karakitsios E, Dokoumetzidis A (2019) PAGE 28 abstr 8895 (Koumaki 2023 reference 9).
  • Description: Two-compartment population PK model for intravenous mecillinam (amdinocillin) in healthy adults, used to determine systemic MIC breakpoints against Enterobacterales for intermittent, extended, and continuous infusion regimens (Koumaki 2023). CL, Q, Vc and Vp are parameterised per kilogram of body weight (linear exponent 1, not allometric). The unbound concentration Cu = fu * Cc (fu = 0.9, i.e. 10 percent protein binding) is the driver of the PK/PD target fT>MIC >= 40 percent of the dosing interval used for probability of target attainment. Parameters were obtained by MCMC (Stan) reanalysis of published AGGREGATE concentration mean and SD data, not individual-level data; no residual error was estimated.
  • Article: Microbiol Spectr 11(2):e03441-22
  • Underlying concentration data: Gambertoglio 1980, Antimicrob Agents Chemother 18:952-956
  • Body-weight reference used by the Monte Carlo simulation: NCHS Series 3 No. 46 (2021)

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 %fT>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).

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%.

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.")
Lognormal body-weight fit vs. the NCHS 2015-2018 percentiles it was derived from.
Sex NCHS (kg)_5th NCHS (kg)_10th NCHS (kg)_15th NCHS (kg)_25th NCHS (kg)_50th NCHS (kg)_75th NCHS (kg)_85th NCHS (kg)_90th NCHS (kg)_95th Fitted lognormal (kg)_5th Fitted lognormal (kg)_10th Fitted lognormal (kg)_15th Fitted lognormal (kg)_25th Fitted lognormal (kg)_50th Fitted lognormal (kg)_75th Fitted lognormal (kg)_85th Fitted lognormal (kg)_90th Fitted lognormal (kg)_95th
Female 49.8 53.9 57.5 62.2 73.1 88.6 98.2 105.3 119.6 48.8 53.7 57.3 63.0 75.2 89.8 98.8 105.4 115.9
Male 61.7 66.6 69.9 75.3 87.4 101.9 110.6 119.4 130.3 61.0 66.2 70.0 75.9 88.5 103.0 111.8 118.2 128.3 

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).")
Fitted lognormal body-weight distributions (NCHS Series 3 No. 46, Tables 3 and 5).
Sex meanlog sdlog Implied mean (kg) NCHS reported mean (kg)
Female 4.3201 0.2631 77.8 77.5
Male 4.4825 0.2261 90.7 90.6

Simulation 

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)
#> ℹ parameter labels from comments will be replaced by 'label()'

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.

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().

# `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.

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%."
)
Simulated vs. published derived PK values (70 kg typical subject). * differs from reference by >20%.
NCA parameter treatment Reference Simulated % diff
CL/F (L/h) 10 mg/kg IV, 15-min infusion 14.5 14.5 -0.0%
Vss (IV) (L) 10 mg/kg IV, 15-min infusion 14.1 14.3 +0.8%

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).

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.")
Derived quantities reported in the text of Koumaki 2023.
Quantity Koumaki 2023 Gambertoglio 1980 (NCA) This model (70 kg, typical)
CL (mL/min/kg) 3.45 3.5 3.45
Vss = V1 + V2 (mL/kg) 202 (text) / 203.68 (V1 + V2) 230 203.68
Terminal half-life (h) not reported not reported 0.739

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% (%fT>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.

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)
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'."
    )
  )
Replicates Table 2 of Koumaki 2023: PTA (%) for intermittent and prolonged infusion. Each cell is ‘published / reproduced’.
MIC (mg/L) 20 min, 1 g x 3 20 min, 1 g x 4 20 min, 1.2 g x 4 2 h, 1 g x 3 2 h, 1 g x 4 2 h, 1.2 g x 4 4 h, 1 g x 3 4 h, 1 g x 4 4 h, 1.2 g x 4
0.125 99.94 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0
0.250 99.66 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0
0.500 98.01 / 100.0 99.98 / 100.0 99.98 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0
1.000 90.39 / 98.0 99.75 / 100.0 99.82 / 100.0 99.97 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0
2.000 63.99 / 74.5 96.55 / 99.5 97.14 / 99.0 98.97 / 99.5 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0
4.000 17.93 / 18.5 72.38 / 82.0 77.09 / 88.5 80.72 / 84.5 99.92 / 100.0 99.94 / 100.0 99.88 / 100.0 100.00 / 100.0 100.00 / 100.0
8.000 0.46 / 0.0 14.37 / 13.5 22.29 / 27.5 12.90 / 14.0 75.27 / 75.0 80.82 / 93.5 64.98 / 67.5 88.26 / 89.5 97.25 / 97.5
16.000 0.00 / 0.0 0.13 / 0.0 0.57 / 0.0 0.05 / 0.0 4.75 / 3.0 11.49 / 13.5 2.20 / 0.5 9.81 / 8.0 25.61 / 24.5

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%.

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)
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'."
    )
  )
Replicates Table 3 of Koumaki 2023: PTA (%) for continuous infusion. Each cell is ‘published / reproduced’.
MIC (mg/L) 150 mg/day 500 mg/day 1,000 mg/day 2,000 mg/day 3,500 mg/day 4,800 mg/day
0.125 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0
0.250 98.53 / 86.5 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0
0.500 33.93 / 5.0 99.91 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0
1.000 0.29 / 0.0 64.81 / 64.5 99.91 / 100.0 100.00 / 100.0 100.00 / 100.0 100.00 / 100.0
2.000 0.00 / 0.0 2.17 / 1.0 64.74 / 63.5 99.91 / 100.0 100.00 / 100.0 100.00 / 100.0
4.000 0.00 / 0.0 0.00 / 0.0 2.16 / 0.5 64.81 / 67.0 99.62 / 100.0 100.00 / 100.0
8.000 0.00 / 0.0 0.00 / 0.0 0.00 / 0.0 2.17 / 0.0 50.17 / 43.5 98.53 / 85.0
16.000 0.00 / 0.0 0.00 / 0.0 0.00 / 0.0 0.00 / 0.0 0.90 / 0.5 33.93 / 7.5

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.

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.")
Continuous-infusion PTA as a function of the dose/MIC ratio.
Dose / MIC ratio Published Reproduced CL cut-off (L/h) Lognormal fit through the ratio 300 and 500 points
218.75 0.90 0.50 8.20 18.3
250.00 2.17 0.50 9.38 24.3
300.00 33.93 6.25 11.25 33.9
437.50 50.17 43.50 16.41 56.8
500.00 64.79 65.00 18.75 64.8
600.00 98.53 85.75 22.50 74.7

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 %fT>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.