---
title: "Meropenem gentamicin ciprofloxacin (Sadouki 2025)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Meropenem gentamicin ciprofloxacin (Sadouki 2025)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

```{r model-meta}
mod_fun <- readModelDb("Sadouki_2025_meropenem_gentamicin_ciprofloxacin")
mod     <- mod_fun()
cat(mod$description, "\n\n", mod$reference, "\n", sep = "")
```

* Article: <https://doi.org/10.1038/s41598-025-29354-y> (Open Access, CC BY 4.0)

This is an **in-vitro pharmacodynamic** model. It is included in
`nlmixr2lib` for reference only; unlike the rest of the library it does not
describe drug exposure in humans or animals. The "subjects" in the source
study are bacterial cultures, and random effects (eta) capture
between-experiment variability rather than between-subject IIV. See
"Assumptions and deviations" at the end of this vignette for the full list of
ways this model departs from the human-pop-PK conventions used elsewhere in
nlmixr2lib.

# Biological context

The source data are from static time-kill experiments against
*Escherichia coli* NCTC 12,241 in 96-well plates (200 uL CAMHB, 37.5 C,
24 h). Eight scaled-MIC levels (0.25, 0.5, 1, 2, 4, 8, 16 x MIC, plus an
antibiotic-free growth control) were tested for each of seven regimens:
three monotherapies (meropenem, gentamicin, ciprofloxacin), three two-way
combinations, and the three-way combination. Two starting inoculum levels
(10^3 and 10^5 CFU/mL) were studied; the lower inoculum reduces the
estimated B0 and Bmax of the logistic growth model.

The experiment-level MICs (used to scale tested concentrations) were
0.03 mg/L for meropenem, 0.015 mg/L for ciprofloxacin, and 1 mg/L for
gentamicin (Sadouki 2025 Methods, Static time-kill experiments).

The full population/experimental metadata is available programmatically:

```{r population}
mod$population
```

# Source trace

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

| Equation / parameter | Value | Source location |
|---|---|---|
| Growth: `dB/dt = (knet*(1 - B/10^Bmax) - sum_drug Effect) * B` | n/a | Methods, Pharmacodynamic modelling (page 2 -- 3) |
| Drug effect: `Effect_drug = Emax * C^h / (IC50^h + C^h) * (1 - beta*(1 - exp(-t*tau)))` | n/a | Methods, Pharmacodynamic modelling (page 3) |
| Mer chemical degradation: `dCmer/dt = -kdeg * Cmer` | n/a | Methods, Pharmacodynamic modelling (page 2) |
| `knet` | 1.35 | Table 1 -- Growth model parameters |
| `B0` | 5.29 log10 CFU/mL | Table 1 |
| `Bmax` | 10 log10 CFU/mL | Table 1 |
| `cat_b0_lowinoc` | -0.326 | Table 1 -- Inoculum effect |
| `cat_bmax_lowinoc` | -0.110 | Table 1 -- Inoculum effect |
| `Emax_Mer`, `IC50_Mer`, `hill_Mer` | 4.18, 0.0781, 2.76 | Table 1 -- Meropenem killing |
| `BETA_Mer`, `TAU_Mer` (FIXED) | 0.922, 0.570 | Table 1 -- Meropenem regrowth |
| `coef_taumer_Mer` | 0.00155 per mg/L Mer | Table 1 -- Proportional MER on TAU |
| `kdeg_Mer` (derived; FIXED) | 0.01317 1/h | Discussion (page 7): "10% drop in 8 h" + Suppl. Fig. 3 |
| `Emax_Gen`, `IC50_Gen`, `hill_Gen` | 5.47, 1.12, 3.63 | Table 1 -- Gentamicin killing |
| `BETA_Gen`, `TAU_Gen` (FIXED) | 0.829, 0.517 | Table 1 -- Gentamicin regrowth |
| `Emax_BetaGen`, `IC50_BetaGen`, `hill_BetaGen` (FIXED) | -2.97, 5.72, 20 | Table 1 -- Emax-on-BETA-Gen |
| `Emax_Cip`, `IC50_Cip`, `hill_Cip` | 4.55, 0.0106, 3.58 | Table 1 -- Ciprofloxacin killing |
| `BETA_Cip`, `TAU_Cip` (FIXED) | 0.674, 0.359 | Table 1 -- Ciprofloxacin regrowth |
| `Emax_BetaCip`, `IC50_BetaCip`, `hill_BetaCip` (FIXED) | -4, 0.017, 20 | Table 1 -- Emax-on-BETA-Cip |
| `combo_beta` (FIXED) | -1 | Table 1 -- Drug interactions, categorical 2/3-way effect |
| `mer_on_ic50cip`, `gen_on_ic50cip` (FIXED) | -0.353, -0.576 | Table 1 -- Drug interactions, proportional effects on Cip IC50 |
| Residual error: additive on log(bacteria) | 0.864 | Table 1 -- Residual variability |

# Virtual experimental grid

The scenarios reproduce the panels of Sadouki 2025 Figure 1: seven
antibiotic-containing regimens at eight scaled MIC levels (0 = growth
control through 16 x MIC). Standard inoculum (10^5 CFU/mL).

```{r grid}
mic <- c(meropenem = 0.03, ciprofloxacin = 0.015, gentamicin = 1)

mic_levels <- c(0, 0.25, 0.5, 1, 2, 4, 8, 16)

regimens <- tibble(
  regimen = c("Growth control",
              "Mer mono", "Gen mono", "Cip mono",
              "Mer + Gen", "Mer + Cip", "Gen + Cip",
              "Mer + Gen + Cip"),
  mer_in  = c(0L, 1L, 0L, 0L, 1L, 1L, 0L, 1L),
  gen_in  = c(0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L),
  cip_in  = c(0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L)
)

scenarios <- tidyr::expand_grid(regimens, mic_mult = mic_levels) |>
  filter(!(regimen == "Growth control" & mic_mult > 0)) |>
  mutate(
    Cmer_init = mer_in * mic_mult * mic[["meropenem"]],
    Cgen      = gen_in * mic_mult * mic[["gentamicin"]],
    Ccip      = cip_in * mic_mult * mic[["ciprofloxacin"]],
    id        = dplyr::row_number()
  )

knitr::kable(scenarios, caption = "Virtual experimental grid (one row per simulation cohort).")
```

# Simulation

Reproduce the published figures with typical-value (eta = 0) trajectories.
The model carries between-experiment etas for visual VPCs but the published
figures are typical-value predictions.

```{r simulate}
mod_typ <- rxode2::zeroRe(mod)

times <- c(seq(0, 8, by = 0.25), seq(9, 24, by = 0.5))

ev <- scenarios |>
  rowwise() |>
  do({
    sc <- .
    tibble(
      id   = sc$id,
      time = c(0, times),
      evid = c(1L, rep(0L, length(times))),
      amt  = c(sc$Cmer_init, rep(NA_real_, length(times))),
      cmt  = c("mer", rep(NA_character_, length(times))),
      CONMED_MER    = sc$mer_in,
      CONMED_GEN    = sc$gen_in,
      CONMED_CIP    = sc$cip_in,
      CONMED_GEN_CC = sc$Cgen,
      CONMED_CIP_CC = sc$Ccip,
      LOWINOC     = 0L,
      regimen     = sc$regimen,
      mic_mult    = sc$mic_mult
    )
  }) |>
  ungroup()

sim <- rxode2::rxSolve(
  mod_typ,
  events     = ev,
  keep       = c("regimen", "mic_mult",
                 "CONMED_MER", "CONMED_GEN", "CONMED_CIP",
                 "CONMED_GEN_CC", "CONMED_CIP_CC"),
  returnType = "data.frame"
) |>
  dplyr::as_tibble() |>
  dplyr::mutate(log10_bact = log10(pmax(bacteria, 1)))
```

# Replicate Figure 1 -- time-kill curves

Each panel of Figure 1 shows log10(CFU/mL) versus time across regimens at
one MIC multiple. The y-axis dashed line at 100 CFU/mL (= 10^2) is the
practical limit of detection.

```{r figure-1, fig.height = 9, fig.width = 9}
fig1_data <- sim |>
  filter(time <= 24) |>
  mutate(panel = paste0(mic_mult, "x MIC"))

fig1_data$panel <- factor(
  fig1_data$panel,
  levels = c("0x MIC", "0.25x MIC", "0.5x MIC", "1x MIC",
             "2x MIC", "4x MIC", "8x MIC", "16x MIC")
)

ggplot(fig1_data, aes(time, log10_bact, colour = regimen)) +
  geom_line(linewidth = 0.7) +
  geom_hline(yintercept = 2, linetype = 2, colour = "grey40") +
  facet_wrap(~ panel, ncol = 3) +
  scale_y_continuous(limits = c(0, 11), breaks = seq(0, 10, by = 2)) +
  labs(
    x = "Time (h)",
    y = "log10 CFU/mL",
    colour = "Regimen",
    title = "Figure 1 -- typical-value time-kill curves",
    caption = "Replicates Figure 1 of Sadouki 2025."
  ) +
  theme_bw() +
  theme(legend.position = "bottom")
```

Sanity checks against Sadouki 2025 Figure 1:

* Growth control panel rises from B0 = 10^5.29 to the carrying capacity
  Bmax = 10^10 by 24 h.
* At 16 x MIC: the three-way (Cip + Gen + Mer) and two-way Cip-containing
  combinations clear bacteria below the limit of detection within 1 -- 2 h
  and stay there; meropenem monotherapy kills early but regrowth emerges by
  24 h (consistent with the unmodified BETA_Mer and the small
  meropenem chemical degradation).
* Gentamicin monotherapy at 8 -- 16 x MIC kills sustainably (no regrowth)
  because the Emax-on-BETA-Gen term suppresses the regrowth shape at high
  Gen concentrations.
* Two-way Mer + Gen at high MIC matches Cip-containing combinations because
  the categorical -1 shift on BETA reverses regrowth in any 2/3-way
  combination.

# Replicate Figure 2 -- concentration -- effect curves

Figure 2 left column: peak killing effect (Effect, 1/h) versus drug
concentration for each drug, separated by mono-therapy and combination
modifiers (combination changes IC50 for ciprofloxacin via the proportional
synergy terms).

```{r figure-2-effect, fig.height = 7, fig.width = 8}
peak_effect <- sim |>
  filter(time == 0) |>  # Effect at t = 0 isolates the Hill curve from the regrowth term
  select(regimen, mic_mult, CONMED_MER, CONMED_GEN, CONMED_CIP,
         eff_mer, eff_gen, eff_cip)

# Hill curves -- peak Emax * C^h / (IC50^h + C^h)
conc_grid <- tibble(
  conc = exp(seq(log(1e-4), log(10), length.out = 200))
)

# Meropenem (constant across regimens for the peak-effect Hill curve)
hill_mer <- conc_grid |>
  mutate(
    eff = 4.18 * conc^2.76 / (0.0781^2.76 + conc^2.76),
    drug = "Meropenem"
  )

# Gentamicin
hill_gen <- conc_grid |>
  mutate(
    eff = 5.47 * conc^3.63 / (1.12^3.63 + conc^3.63),
    drug = "Gentamicin"
  )

# Ciprofloxacin -- one curve per IC50 multiplier (mono / +Mer / +Gen / +Mer+Gen)
ic50_cip_base <- 0.0106
ic50_cip_mods <- tibble(
  presence = c("Cip alone", "Cip + Mer", "Cip + Gen", "Cip + Mer + Gen"),
  factor   = c(1,
               1 + (-0.353),
               1 + (-0.576),
               1 + (-0.353) + (-0.576))
)
hill_cip <- tidyr::expand_grid(conc_grid, ic50_cip_mods) |>
  mutate(
    ic50_eff = ic50_cip_base * factor,
    eff      = 4.55 * conc^3.58 / (ic50_eff^3.58 + conc^3.58),
    drug     = paste0("Ciprofloxacin (", presence, ")")
  )

ggplot() +
  geom_line(data = hill_mer, aes(conc, eff), colour = "steelblue", linewidth = 0.8) +
  geom_line(data = hill_gen, aes(conc, eff), colour = "darkgreen", linewidth = 0.8) +
  geom_line(data = hill_cip, aes(conc, eff, colour = presence), linewidth = 0.7) +
  facet_wrap(~ NULL) +
  scale_x_log10() +
  labs(
    x = "Drug concentration (mg/L, log scale)",
    y = "Peak killing effect (1/h)",
    colour = "Cip context",
    title = "Figure 2 (left column) -- concentration -- effect curves",
    caption = paste("Replicates Figure 2 of Sadouki 2025;",
                    "blue = Mer, green = Gen, multi-coloured = Cip variants.")
  ) +
  theme_bw() +
  theme(legend.position = "bottom")
```

Note the leftward shift of the ciprofloxacin curve when meropenem and / or
gentamicin are present -- this is the synergy on potency (lower IC50_Cip)
quantified by `mer_on_ic50cip = -0.353` and `gen_on_ic50cip = -0.576` in
Table 1.

# Replicate Figure 2 -- loss-of-susceptibility curves

Figure 2 right column: relative effect over time (Effect(t) / Emax) at
several scaled concentrations. The model's regrowth term
`(1 - beta * (1 - exp(-t * tau)))` decays from 1 (full effect at t = 0) to
`(1 - beta)` at large t. When `beta < 0` (combination case, or high
Gen / Cip concentrations triggering the Emax-on-BETA term), the term grows
and effect is *enhanced* over time.

```{r figure-2-loss, fig.height = 6, fig.width = 8}
t_grid <- seq(0, 24, by = 0.1)

loss_curve <- function(beta, tau) {
  1 - beta * (1 - exp(-t_grid * tau))
}

# Meropenem panel: monotherapy (BETA = 0.922) and combination (BETA = -0.078)
mer_curves <- bind_rows(
  tibble(time = t_grid, rel_eff = loss_curve(0.922, 0.570),
         drug = "Mer", context = "monotherapy"),
  tibble(time = t_grid, rel_eff = loss_curve(0.922 - 1, 0.570),
         drug = "Mer", context = "in 2/3-way combination")
)

# Gentamicin: typical mono BETA = 0.829; concentration-dependent BETA shift
# at 0.25, 1, 8, 16 x MIC; combination case
gen_concs <- c(0.25, 1, 8, 16) * mic[["gentamicin"]]
beta_gen_shift <- function(c, ic50 = 5.72, emax = -2.97, h = 20) {
  emax * c^h / (ic50^h + c^h)
}
gen_curves <- bind_rows(
  lapply(gen_concs, function(c) {
    tibble(time = t_grid,
           rel_eff = loss_curve(0.829 + beta_gen_shift(c), 0.517),
           drug = "Gen",
           context = paste0(c, " mg/L (mono)"))
  })
) |>
  bind_rows(
    tibble(time = t_grid,
           rel_eff = loss_curve(0.829 - 1, 0.517),
           drug = "Gen",
           context = "in 2/3-way combination")
  )

# Ciprofloxacin: similar at several concentrations
cip_concs <- c(0.25, 1, 4, 16) * mic[["ciprofloxacin"]]
beta_cip_shift <- function(c, ic50 = 0.017, emax = -4, h = 20) {
  emax * c^h / (ic50^h + c^h)
}
cip_curves <- bind_rows(
  lapply(cip_concs, function(c) {
    tibble(time = t_grid,
           rel_eff = loss_curve(0.674 + beta_cip_shift(c), 0.359),
           drug = "Cip",
           context = paste0(c, " mg/L (mono)"))
  })
) |>
  bind_rows(
    tibble(time = t_grid,
           rel_eff = loss_curve(0.674 - 1, 0.359),
           drug = "Cip",
           context = "in 2/3-way combination")
  )

loss_data <- bind_rows(mer_curves, gen_curves, cip_curves) |>
  mutate(drug = factor(drug, levels = c("Mer", "Gen", "Cip")))

ggplot(loss_data, aes(time, rel_eff, colour = context)) +
  geom_line(linewidth = 0.7) +
  geom_hline(yintercept = 1, linetype = 3, colour = "grey50") +
  facet_wrap(~ drug, ncol = 1, scales = "free_y") +
  labs(
    x = "Time (h)", y = "Relative effect (1 - BETA*(1-exp(-t*TAU)))",
    colour = "Context",
    title = "Figure 2 (right column) -- loss / enhancement of effect over time",
    caption = "Replicates Figure 2 (right) of Sadouki 2025."
  ) +
  theme_bw() +
  theme(legend.position = "right")
```

# Validation against published findings

The model file faithfully reproduces the qualitative behaviours called out
in Sadouki 2025:

* **Sustained killing in any 2/3-way combination** -- the categorical -1
  shift on BETA reverses regrowth so that effect is sustained rather than
  waning. Visible in the 1 -- 16 x MIC panels of the simulated Figure 1.
* **Synergy on potency for ciprofloxacin** -- the simulated peak-effect
  Hill curves for Cip shift left when meropenem and / or gentamicin are
  present (Figure 2 left). The 35.3% IC50 reduction (MER on Cip) and
  57.6% reduction (GEN on Cip) match Table 1.
* **Three-way vs two-way equivalence** -- in Figure 1, three-way
  Mer + Gen + Cip and two-way Cip + Mer / Cip + Gen at high MIC have
  near-identical kill trajectories, because the dominant interaction terms
  on BETA and IC50_Cip already saturate at the two-way level.
* **Mer + Gen indifference on cidal effects** -- at the per-drug Hill-curve
  level the Mer + Gen combination has no proportional effect on Mer or Gen
  IC50; the only modifier is the categorical -1 BETA shift that prevents
  regrowth (Sadouki 2025 Discussion, page 6).

A side-by-side numerical NCA comparison (Cmax, AUC, half-life) is not
relevant here -- there is no PK exposure profile to integrate. PKNCA is
intentionally **not** used in this vignette; see "Assumptions and
deviations" below.

# Assumptions and deviations

This model deviates from several `nlmixr2lib` conventions because the
underlying paper is an **in-vitro bacterial pharmacodynamic** study, not a
human / animal population PK / PD study. The deviations were sanctioned by
the operator at sidecar time (task `032-husain_135`, request 001
response: option B, "extract anyway with clear in-vitro metadata, file
under `inst/modeldb/pharmacodynamics/`").

* **Compartments are not the canonical `depot` / `central` / `effect` / `Cc`
  set.** The state variables are `mer` (meropenem solution concentration in
  mg/L) and `bacteria` (E. coli concentration in CFU/mL on linear scale).
  Bacterial counts are observed as `log_bact = log(bacteria)`. Gentamicin
  and ciprofloxacin concentrations are passed as static covariates (`Cgen`,
  `Ccip`) because the experimental design holds them constant for 24 h.
* **Covariate columns are not in the canonical
  `inst/references/covariate-columns.md` register.** The register applies
  to human pop-PK covariates (body weight, renal function, etc.). Here the
  covariates are experimental-design indicators (`MER_PRESENT`,
  `GEN_PRESENT`, `CIP_PRESENT`, `LOWINOC`) and static drug concentrations
  (`Cgen`, `Ccip`). Each is documented in the model file's `covariateData`
  list with its `source_name`, but no canonical-register entry was created.
* **Population metadata uses in-vitro fields.** The standard fields
  (`age_range`, `weight_range`, `sex_female_pct`, etc.) are not applicable.
  The model file's `population` list instead records `organism`, `medium`,
  `temperature`, `duration`, `mic_values`, `concentration_range`,
  `inoculum_options`, and `regimens`.
* **Random effects are between-experiment, not between-subject.** The paper's
  Methods section (page 2): "Random effects (eta) were distinguishing
  variability between experiments". The eta values in `ini()` are kept for
  faithfulness to the published model but are not "IIV" in the standard
  pop-PK sense.
* **PKNCA validation is replaced with figure replication.** PKNCA is not
  the right validation target for an in-vitro bacterial-load model -- there
  is no dose, no absorption, no exposure profile to integrate. The
  vignette instead reproduces Sadouki 2025 Figure 1 (kill curves over 24 h
  across all seven regimens and eight MIC multiples) and Figure 2 (Hill
  concentration -- effect curves and loss-of-susceptibility curves).
* **Parameter-naming exceptions to `naming-conventions.md`.** Drug-specific
  parameters use a `<param>_<drug>` suffix pattern (e.g., `lemax_mer`,
  `ic50_cip_eff`) rather than the canonical `lcl` / `lvc` PK-parameter
  names. Parameter `etabeta_mer` is on the linear scale (matching the
  paper's "additive eta on BETA" footnote in Table 1). These deviations
  match the structure of the published model and are mandatory for a
  faithful reproduction.
* **Meropenem chemical degradation rate `kdeg_Mer = 0.01317 1/h` is
  derived, not directly tabulated.** The paper's Discussion (page 7) and
  Suppl. Fig. 3 report "10% drop in 8 h"; this skill computes
  `kdeg = -log(0.9) / 8 = 0.01317 1/h` and fixes the value. Recorded as a
  derived value with provenance in the model file's in-line comment.
* **The supplement (Suppl. Fig. 1 -- 3, Suppl. Tables) was not on disk** at
  extraction time. Parameter values are taken from Table 1 of the main
  paper, which the paper explicitly designates as the "Final parameter
  estimates of the pharmacodynamic model". No initial estimates were used.

# Errata

No errata or corrigenda are listed on the journal landing page or in
PubMed search as of 2026-05-08.
