---
title: "Meropenem (Bergen 2017)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Meropenem (Bergen 2017)}
  %\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

* Citation: `r rxode2::rxode(readModelDb("Bergen_2017_meropenem"))$reference`
* Description: `r rxode2::rxode(readModelDb("Bergen_2017_meropenem"))$description`
* Article: <https://doi.org/10.1128/AAC.02642-16>

This is **not** a population PK model. It is a mechanism-based PK/PD model (MBM)
of bacterial killing and resistance, fit with S-ADAPT to viable counts for
*Pseudomonas aeruginosa* 1280 in the dynamic hollow-fiber infection model
(HFIM) across three simulated renal-function scenarios (augmented renal
clearance, normal, impaired) and four dosing regimens. Meropenem exposure is
reproduced as a concentration state `cmem` that the user doses and that
declines with a fixed, renal-function-specific half-life (the underlying
disposition came from a 2-compartment popPK model in critically ill patients
that the paper cites as reference 20 / Mattioli 2016 but does not republish in
full). Because there is no clinical absorption-distribution-elimination
profile of a single drug to integrate for an NCA, NCA / PKNCA is not an
appropriate validation; the checks below are the mechanistic equivalents
(carrying-capacity hold, mass-balance behaviour, replicate of the published
Figure 2 kill / regrowth trajectories, and a comparison against Table 4 PK
summary endpoints).

# Population (biological context)

The model describes *P. aeruginosa* 1280, a meropenem-susceptible clinical
isolate (Etest MIC 0.25 mg/L) from a critically ill patient with a
soft-tissue infection. It was grown in cation-adjusted Mueller-Hinton broth
in the HFIM at 36 C for 10 days. Three renal-function scenarios were
simulated, with meropenem exposure profiles generated from the Mattioli 2016
critically-ill popPK (paper reference 20):

| Scenario | CLcr (mL/min) | Meropenem CL (L/h) | t1/2 (h) |
|---|---|---|---|
| Augmented renal clearance (ARC) | ~285 | 34.0 | 0.6 |
| Normal renal function           | 120  | 16.3 | 1.1 |
| Impaired renal function         | ~10  | 4.1  | 4.0 |

Four dosing regimens were studied (2 g q8h, 1 g q8h, 0.5 g q8h, all as
30-min IV infusions; plus 1 g q12h administered to the impaired-renal-function
scenario only). All MBM parameters in the model file are the typical-value
estimates from Bergen 2017 Table 3; the simulated renal-function half-life
defaults to 1.1 h (normal renal function) and can be overridden at
simulation time.

The same information is available programmatically via
`readModelDb("Bergen_2017_meropenem")$population`.

# Source trace

Per-parameter origins are recorded as in-file comments next to each `ini()`
entry in `inst/modeldb/specificDrugs/Bergen_2017_meropenem.R`. All parameter
values are the HFIM population-mean estimates from Bergen 2017 Table 3; the
PK summaries (Cmax, Cmin, AUC, %fT>MIC) are from Bergen 2017 Table 4.

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `log10cfu0` (initial inoculum) | 6.97 log10 CFU/mL | Table 3 (Log CFU0) |
| `log10cfumax` (max population size) | 9.98 log10 CFU/mL | Table 3 (Log CFUmax) |
| `lk21` (replication rate, FIXED) | 50 /h | Table 3 (footnote a) |
| `mgt_s`, `mgt_i`, `mgt_r` (mean generation times) | 49.3, 683, 78.5 min | Table 3 (k12,X^-1 rows; k12 = 60/MGT) |
| `log10mf_i`, `log10mf_r` (mutation freqs) | -3.66, -6.28 | Table 3 (LogMF_I, LogMF_R) |
| `imax_rep_s` (FIXED), `imax_rep_i`, `imax_rep_r` | 1.0, 0.673, 0.956 | Table 3 (Imax_Rep_X; footnote b for S) |
| `ic50_rep_s`, `ic50_rep_i`, `ic50_rep_r` | 0.648, 2.96, 6.09 mg/L | Table 3 (IC50_Rep_X) |
| `hill_s` (FIXED), `hill_i`, `hill_r` | 1.0, 7.14, 2.19 | Table 3 (Hill_X; footnote c for S) |
| `addSd` (residual SD, log10 scale) | 0.493 | Table 3 (SD_CFU) |
| `thalf_mem` (FIXED) | 1.1 h (default; 0.6/4.0 for ARC/impaired) | Table 4 (t_1/2 per renal function) |
| Life-cycle growth ODEs (2 states/subpop) | n/a | Methods (Eqs around 5-6 in the prose) |
| Replication factor `REP = 2 * (1 - CFUall/CFUmax)` | n/a | Methods (Eq 3) |
| Inhibition of replication (Hill) | n/a | Methods (Eq 4) |

## Units (dimensional analysis)

| Symbol | Meaning | Units |
|--------|---------|-------|
| `bact_susceptible1`, `bact_susceptible2`, `bact_intermediate1`, `bact_intermediate2`, `bact_resistant1`, `bact_resistant2` | bacterial states | CFU/mL |
| `cmem` | meropenem concentration | mg/L |
| `lk21`, `k12*`, `kel_mem` | rate constants | 1/h |
| `ic50_rep_*` | half-effect concentrations | mg/L |
| `hill_*`, `rep_factor`, `irep_*`, `imax_rep_*` | dimensionless | -- |
| `cfumax`, `cfu0`, `cfu_all`, `cfu_less_susc` | population scale / inoculum / sums | CFU/mL |

Every growth ODE term has the form (1/h) * (CFU/mL) = (CFU/mL)/h, matching
`d/dt(state)`; the meropenem ODE has (1/h) * (mg/L) = (mg/L)/h. `k12 = 60/MGT`
carries a hidden 60 in (min/h), converting the mean generation time (min) to
a rate (1/h).

```{r load-model}
mod <- rxode2::rxode(readModelDb("Bergen_2017_meropenem"))
mod$state
```

# Parameter table (paper vs. file)

```{r param-table}
params <- mod$theta
knitr::kable(
  data.frame(parameter = names(params), file_value = unname(params)),
  caption = "Typical-value parameters loaded from the model file (Bergen 2017 Table 3 + Table 4 default)."
)
```

# Solver settings used in this vignette

The two-state life-cycle model is moderately stiff because `k21` (50 /h) is
much faster than the subpopulation growth rate constants `k12` (~0.09-1.2 /h).
Default lsoda tolerances accumulate enough error over ~24 h of growth to drive
the trajectory to `NaN`. Tighten the absolute and relative tolerances when
solving:

```{r solver-opts}
solver_opts <- list(maxsteps = 1e6, atol = 1e-12, rtol = 1e-8)
```

# Carrying-capacity (growth control) check

With no antibiotic, the population grows from the inoculum (10^6.97) toward
a stationary plateau. As in Rees 2018 (the cousin two-state MBM by the same
group), the published equations have a structural property:
`rep_factor = 2 * (1 - cfu_all / cfumax)` settles at **1** at steady state
(each surviving division replaces one cell). Solving for `rep_factor = 1`
gives `cfu_all = 0.5 * cfumax`, so the realized stationary population is
`0.5 * 10^9.98`, i.e. `log10(0.5 * 10^9.98) = 9.679` log10 CFU/mL -- not 9.98.
This matches the growth-control plateau visible in Bergen 2017 Figure 2
(~9.5-10 log10 CFU/mL) and is correct behaviour of the published equations
rather than a transcription error.

```{r growth-control}
ev_gc <- as.data.frame(rxode2::et(time = seq(0, 96, by = 0.5)))
gc <- do.call(rxode2::rxSolve, c(list(object = mod, events = ev_gc,
                                      returnType = "data.frame"), solver_opts))

cat(sprintf("Inoculum log10CFU(0)  = %.3f  (Table 3 Log10CFU0 = 6.97)\n", gc$Cc[1]))
cat(sprintf("Plateau  log10CFU(96) = %.3f  (expected 0.5*CFUmax = %.3f)\n",
            tail(gc$Cc, 1), log10(0.5) + 9.98))

ggplot(gc, aes(time, Cc)) +
  geom_line(linewidth = 1) +
  geom_hline(yintercept = log10(0.5) + 9.98, linetype = 2, colour = "grey50") +
  labs(x = "Time (h)", y = expression(log[10] ~ CFU/mL),
       title = "Antibiotic-free growth control",
       caption = "Dashed: stationary population 0.5 * CFUmax = 9.679 log10 CFU/mL.")
```

# Steady-state IV-infusion dose helper

The model doses the meropenem **concentration** `cmem` (mg/L) directly, the
same idiom Rees 2018 uses. A 0.5-h IV infusion every tau hours that reaches a
prescribed steady-state Cmax has rate `R = Cmax * kel * (1 - exp(-kel*tau)) /
(1 - exp(-kel*Tinf))` and per-dose amount `R * Tinf`. We use the published
Table 4 unbound Cmax to set the rate for each regimen so the simulated `cmem`
profile reproduces the paper's exposure summaries.

```{r rate-helper}
rate_for_cmax <- function(Cmax, kel, tau, Tinf) {
  Cmax * kel * (1 - exp(-kel * tau)) / (1 - exp(-kel * Tinf))
}
```

# Replicate Figure 2 (HFIM kill / regrowth trajectories)

Bergen 2017 Figure 2 shows total viable counts over 10 days for three
dosing regimens (2, 1, 0.5 g q8h) at each of the three renal functions, plus
a 1 g q12h regimen for the impaired-renal-function scenario, alongside a
growth control. We reproduce the qualitative behaviour reported in the
Results / Table 4 "Outcome" column:

* **ARC (CLcr ~285 mL/min):** every regimen regrows ("RR") -- minimal initial
  killing followed by replacement by less-susceptible bacteria.
* **Normal renal function (CLcr 120 mL/min):** 0.5 g and 1 g q8h regrow
  ("RR"); 2 g q8h suppresses regrowth ("SR").
* **Impaired renal function (CLcr ~10 mL/min):** all regimens (including
  1 g q12h) suppress regrowth ("SR").

```{r fig2-scenarios}
# Table 4 unbound steady-state Cmax (mg/L) for each (renal function, regimen).
# Used to set the infusion rate so cmem mirrors the published exposures.
scenarios <- tribble(
  ~scenario,                  ~thalf,  ~regimen,        ~dose_g, ~Cmax,  ~tau, ~outcome,
  "ARC (CLcr 285)",           0.6,     "2 g q8h",       2.0,     51.3,   8,    "RR",
  "ARC (CLcr 285)",           0.6,     "1 g q8h",       1.0,     25.7,   8,    "RR",
  "ARC (CLcr 285)",           0.6,     "0.5 g q8h",     0.5,     12.8,   8,    "RR",
  "Normal (CLcr 120)",        1.1,     "2 g q8h",       2.0,     70.1,   8,    "SR",
  "Normal (CLcr 120)",        1.1,     "1 g q8h",       1.0,     35.0,   8,    "RR",
  "Normal (CLcr 120)",        1.1,     "0.5 g q8h",     0.5,     17.5,   8,    "RR",
  "Impaired (CLcr ~10)",      4.0,     "2 g q8h",       2.0,     114.5,  8,    "SR",
  "Impaired (CLcr ~10)",      4.0,     "1 g q8h",       1.0,     57.2,   8,    "SR",
  "Impaired (CLcr ~10)",      4.0,     "0.5 g q8h",     0.5,     28.6,   8,    "SR",
  "Impaired (CLcr ~10)",      4.0,     "1 g q12h",      1.0,     49.1,  12,    "SR"
)
scenarios$renal <- factor(scenarios$scenario,
                          levels = c("ARC (CLcr 285)", "Normal (CLcr 120)",
                                     "Impaired (CLcr ~10)"))
```

```{r fig2-simulate}
Tinf <- 0.5  # infusion duration (h)
TMAX <- 240  # 10 days of follow-up

simulate_scenario <- function(thalf, tau, Cmax, id_offset) {
  kel <- log(2) / thalf
  R   <- rate_for_cmax(Cmax, kel, tau, Tinf)
  dose_ev <- as.data.frame(
    rxode2::et(id = id_offset, amt = R * Tinf, rate = R,
               ii = tau, cmt = "cmem", until = TMAX)
  )
  obs_ev <- as.data.frame(
    rxode2::et(id = id_offset, time = seq(0, TMAX, by = 1))
  )
  ev <- bind_rows(dose_ev, obs_ev) %>% arrange(id, time)
  out <- do.call(rxode2::rxSolve,
                 c(list(object = mod, events = ev,
                        params = c(thalf_mem = thalf),
                        returnType = "data.frame"), solver_opts))
  out
}

# Growth controls (one per renal function): no meropenem dosing
gc_runs <- lapply(seq_len(3), function(i) {
  thalf <- c(0.6, 1.1, 4.0)[i]
  scn   <- levels(scenarios$renal)[i]
  obs <- as.data.frame(rxode2::et(id = 1000 + i, time = seq(0, TMAX, by = 2)))
  out <- do.call(rxode2::rxSolve,
                 c(list(object = mod, events = obs,
                        params = c(thalf_mem = thalf),
                        returnType = "data.frame"), solver_opts))
  out$scenario <- scn
  out$regimen  <- "Growth control"
  out
})

trt_runs <- lapply(seq_len(nrow(scenarios)), function(i) {
  r <- scenarios[i, ]
  out <- simulate_scenario(thalf = r$thalf, tau = r$tau, Cmax = r$Cmax,
                           id_offset = i)
  out$scenario <- as.character(r$scenario)
  out$regimen  <- r$regimen
  out
})

sim <- bind_rows(c(gc_runs, trt_runs))
sim$renal <- factor(sim$scenario,
                    levels = c("ARC (CLcr 285)", "Normal (CLcr 120)",
                               "Impaired (CLcr ~10)"))
sim$regimen <- factor(sim$regimen,
                      levels = c("Growth control", "2 g q8h", "1 g q8h",
                                 "0.5 g q8h", "1 g q12h"))
```

```{r fig2-plot, fig.height = 6.5}
loc_total <- 1.0  # limit of counting on antibiotic-free agar (paper Methods)

plot_df <- sim %>%
  mutate(log10_total = pmax(log10(pmax(cfu_all, 1e-6)), loc_total))

ggplot(plot_df, aes(time, log10_total, colour = regimen)) +
  geom_hline(yintercept = loc_total, linetype = 3, colour = "grey60") +
  geom_line(linewidth = 0.7) +
  facet_wrap(~renal, ncol = 1) +
  scale_x_continuous(breaks = seq(0, TMAX, by = 48)) +
  scale_y_continuous(limits = c(0, 10.5)) +
  labs(x = "Time (h)", y = expression(log[10] ~ CFU/mL ~ "(total viable count)"),
       colour = NULL,
       title = "Replicates Bergen 2017 Figure 2",
       caption = paste("Dotted: limit of counting (1.0 log10 CFU/mL).",
                       "Growth control plateau ~9.68 log10 CFU/mL = 0.5 x CFUmax.")) +
  theme(legend.position = "bottom")
```

```{r fig2-endpoints}
sim %>%
  group_by(renal, regimen) %>%
  summarise(
    nadir_total = round(min(log10(pmax(cfu_all, 1e-6))), 2),
    end_total   = round(log10(pmax(tail(cfu_all, 1), 1e-6)), 2),
    end_less_susc = round(log10(pmax(tail(cfu_less_susc, 1), 1e-6)), 2),
    .groups = "drop"
  ) %>%
  knitr::kable(
    caption = "Simulated nadir and day-10 (240 h) total + less-susceptible (intermediate + resistant) populations by regimen and renal function. ARC regimens regrow to the ~9.68 log10 plateau; normal-renal 2 g q8h and all impaired-renal regimens hold the bacteria near the limit of counting.")
```

# PK exposure comparison against Bergen 2017 Table 4

The simulated `cmem` trajectory should reproduce the published Cmax / Cmin
for each scenario (because the dosing was set by inverting the
steady-state IV-infusion formula). We extract the simulated Cmax (end of
infusion) and Cmin (immediately before next dose) for each scenario from
the simulation above and put them next to the Table 4 values.

```{r pk-summary}
# For each (scenario, regimen), find Cmax (max cmem in last full interval)
# and Cmin (cmem at end of last full interval).
pk_obs <- sim %>%
  filter(regimen != "Growth control") %>%
  group_by(renal, regimen) %>%
  filter(time >= TMAX - 24) %>%   # last day at steady state
  summarise(Cmax_sim = max(cmem), Cmin_sim = min(cmem), .groups = "drop")

pk_paper <- scenarios %>%
  transmute(renal, regimen = factor(regimen,
                                    levels = c("2 g q8h", "1 g q8h",
                                               "0.5 g q8h", "1 g q12h")),
            Cmax_paper = Cmax,
            # Table 4 reported Cmin values
            Cmin_paper = c(0.01, 0.00, 0.00,
                           0.50, 0.25, 0.12,
                           31.5, 15.7, 7.9, 6.8))

cmp <- left_join(pk_paper, pk_obs, by = c("renal", "regimen")) %>%
  mutate(Cmax_pct_err = round(100 * (Cmax_sim - Cmax_paper) / Cmax_paper, 1))

knitr::kable(cmp, digits = 2,
             caption = "Simulated steady-state Cmax / Cmin (mg/L) vs Bergen 2017 Table 4. Cmax should match within ~1% by construction (the dose was inverted from the target); Cmin agreement depends on the 1-compartment approximation matching the upstream 2-compartment popPK well enough for the simulated half-life.")
```

# Assumptions and deviations

* **Model class / species.** This is an in-vitro mechanism-based PK/PD model,
  not a popPK model; `population$species` records *P. aeruginosa* 1280. No
  PKNCA validation is performed -- there is no clinical PK profile to
  integrate. The mechanistic checks above (carrying-capacity hold, Figure 2
  replication, Table 4 exposure comparison) replace it, per the
  endogenous/mechanistic validation strategy in the extraction skill.
* **File naming.** The dispatch metadata listed the drug as "Antimicrobial
  Agents and Chemo", which is the journal name (Antimicrobial Agents and
  Chemotherapy), not a drug. The paper unambiguously models meropenem, so
  the model file and this vignette use `Bergen_2017_meropenem`.
* **Meropenem disposition.** The simulated meropenem half-life in the HFIM
  (0.6 / 1.1 / 4.0 h for ARC / normal / impaired renal function) is a fixed
  input taken from the upstream Mattioli 2016 popPK in critically ill
  patients (paper reference 20) and reported in Table 4. It is not an MBM
  estimate. The default in the model file is 1.1 h (normal renal function);
  the vignette overrides it for ARC and impaired-renal-function simulations.
* **One-compartment approximation.** The upstream popPK is a two-compartment
  model in critically ill patients; here we approximate the simulated HFIM
  concentration-time profile with a single first-order elimination state
  parameterised by `thalf_mem`. The Cmax match is exact by construction (the
  per-regimen rate is back-computed from the Table 4 target Cmax) but the
  Cmin agreement depends on how close the 1-compartment exponential decline
  is to the upstream 2-compartment profile, which can differ when there is
  meaningful distributional clearance. Use this model file for MBM behaviour;
  use a 2-compartment PK model if you need exact Cmin reproduction.
* **Mean generation time.** Table 3 labels the growth rows "k12,X^-1 (min)";
  these are mean generation times in minutes, not rate constants. The
  supplement defines `k12 = 60/MGT` (1/h). The model file stores the MGTs
  and derives `k12` in `model()`.
* **Stationary population is 0.5 x CFUmax.** A structural property of the
  published two-state model: at steady state `rep_factor = 1` requires
  `cfu_all = 0.5 * cfumax`, so the realized plateau is at
  log10(0.5 * 10^9.98) = 9.679, not 9.98. The growth-control trajectory above
  hits 9.679 within 24 h and matches the ~9.5-10 log10 CFU/mL plateau in
  Bergen 2017 Figure 2.
* **Mechanism is replication inhibition only.** Bergen 2017 explicitly
  reports that "A model incorporating inhibition of successful replication
  of all three bacterial subpopulations by meropenem best described the
  antibacterial effect. No additional functions or model complexities were
  necessary." There is no direct-killing term in the published model and
  none is included here, in contrast to the Rees 2018 cousin model which
  uses a Hill direct-killing term.
* **Subpopulations are mechanistic, not plate-defined.** The MBM's CFU_S,
  CFU_I, CFU_R are estimated subpopulations seeded by `log10mf_i` /
  `log10mf_r`; the Discussion states "the estimated subpopulations did not
  directly reflect bacterial counts on meropenem-containing agar plates at
  5x and 10x MIC." The derived `cfu_less_susc` output is intermediate +
  resistant, but should not be used as a direct surrogate for agar-plate
  resistance counts.
* **Initial-state partitioning.** All bacteria start in life-cycle state 1;
  state-2 initial conditions are 0 (paper Methods). The balanced-growth
  state-2 fraction is ~k12/k21 (~2-3% for susceptible, < 1% for intermediate
  and resistant), so this introduces only a brief start-up transient that
  matters for the first 5-10 min but is irrelevant on the 10-day time scale.
* **Solver tolerances.** Default lsoda settings (atol=1e-8, rtol=1e-6) drive
  the trajectory to NaN within ~3 h because the population spans 10^6 ->
  10^10 CFU/mL and the relative-error budget is overwhelmed. The vignette
  uses `atol = 1e-12, rtol = 1e-8, maxsteps = 1e6`; users should pass the
  same tolerances to `rxSolve`.
* **Below limit of counting.** Displayed total counts are floored at the
  published limit of counting (1.0 log10 CFU/mL); the underlying states are
  not floored. The limit on antibiotic-containing agar is 0.7 log10 CFU/mL
  per the paper Methods (200 uL of sample plated to increase sensitivity);
  the less-susceptible derived output above is not floored because it does
  not directly correspond to an agar-plate count.
* **No IIV, no covariates.** Bergen 2017 is an in vitro experiment with a
  single isolate. The model has no eta parameters and no covariate columns
  -- `covariateData = list()`. Between-curve variability of 15% CV was set
  during S-ADAPT estimation (Methods); this is an estimation control rather
  than a deployable random effect and is not encoded in the model file.
* **Convention deviations** (`checkModelConventions()` warnings, no
  errors). All are expected for an in-vitro mechanism-based model:
  (a) the bacterial-state compartments (`bact_susceptible1`/`bact_susceptible2`,
  `bact_intermediate1`/`bact_intermediate2`, `bact_resistant1`/`bact_resistant2`)
  and antibiotic-concentration state (`cmem`) are mechanism-specific;
  (b) `lk21` is a fixed mechanistic rate constant that is log-transformed
  for parameterisation consistency; (c) the single observation `Cc` carries
  a non-PK output (log10 viable count, not a drug concentration); (d) the
  dosing/concentration units are `mg/L` because the antibiotic input is a
  concentration in the in-vitro system.
