---
title: "Gentamicin bactericidal activity and adaptive resistance (Mohamed 2012)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Gentamicin bactericidal activity and adaptive resistance (Mohamed 2012)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

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

This is **not** a population PK model. It is a semi-mechanistic PKPD model of
gentamicin bactericidal activity against *Escherichia coli* ATCC 25922, fit in
NONMEM (Laplacian, ADVAN9) to bacterial-count data (natural log of CFU/mL)
from 48 static and 25 dynamic in-vitro time-kill experiments (Mohamed 2012
Methods, Results). The gentamicin exposure is reproduced here as a single
concentration state (`cgent`, mg/L) that the user doses; the structural
neonatal popPK described in the paper's Discussion (and used for the
clinical-dosing predictions in Figs 4-6) is from a separate upstream paper
(Nielsen 2011 / ref 35) and is not packaged here. Because there is no
absorption-distribution-elimination PK profile to integrate for an NCA,
PKNCA is not an appropriate validation; the checks below are the mechanistic
equivalents (carrying-capacity hold, static and dynamic kill replication, and
the adaptive-resistance time-course arithmetic).

# Population (biological context)

The model describes *E. coli* ATCC 25922 (gentamicin MIC 2 mg/L by macrodilution)
grown in cation-adjusted Mueller-Hinton broth at 35 C. Two experimental designs
were used (Mohamed 2012 Methods): **static** time-kill curves at constant
gentamicin concentrations of 0.125, 0.25, 0.5, 1, 2, 4, and 16 mg/L (started at
inoculum ~5e5 CFU/mL, plus a high-inoculum subset pre-grown 12 h to ~1e9 CFU/mL
and exposed to 1, 2, 4, or 16 mg/L), and **dynamic** time-kill curves in a
two-compartment in-vitro kinetic flask simulating preterm-neonate PK with peak
concentrations of 2.0, 3.9, 7.8, and 16 mg/L (matching neonatal 1, 2, 4, and
8 mg/kg doses) as single or repeated doses. The simulated neonatal kinetics use
two exchanger flow rates corresponding to gentamicin rate constants of 0.33 /h
(first 4 h) and 0.037 /h thereafter.

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

# Source trace

Per-parameter origins are recorded as in-file comments next to each `ini()`
entry in `inst/modeldb/specificDrugs/Mohamed_2012_gentamicin.R`. All numeric
values are the final population-mean estimates from Mohamed 2012 Table 1.
The model structure follows Figure 1 and Equations 1-7 of the main text.

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `kgrowth` (bacterial growth rate constant) | 2.00 /h | Table 1 |
| `kdeath` (natural death rate constant; fixed) | 0.179 /h | Table 1 (fix; from Nielsen 2007 ref 36) |
| `bp` (breakpoint for non-zero S->R transfer) | 2.09e6 CFU/mL | Table 1 |
| `bmax` (stationary-phase total count) | 8.26e8 CFU/mL | Table 1 |
| `emax0` (max kill rate at zero AR) | 51.0 /h | Table 1 |
| `ec50` (gentamicin EC50) | 9.93 mg/L | Table 1 |
| `ar50` (AR_on at half-max E_max inhibition) | 0.113 | Table 1 |
| `kon` (AR development rate constant) | 0.0426 L/(mg*h) | Table 1 |
| `koff` (return-to-susceptibility rate constant; fixed) | 0.0139 /h | Table 1 (fix) |
| `s0` (initial inoculum used in predictions) | 4.83e5 CFU/mL | Methods (average across experiments) |
| `lkel` (log in-vitro flask elimination; fixed input) | log(0.037) default | Methods (set to 0 for static, 0.33 then 0.037 /h for dynamic) |
| `addSd` (residual SD on ln(CFU/mL)) | 1.69 | Table 1 (RE_static); RE_dynamic = 2.80 and RRE = 0.618 also reported |
| ODEs for bact_growing / bact_resting (Eqs 1 + 5) | n/a | Main text Eqs 1, 5 |
| Adaptive-resistance binding ODEs (Eqs 6-7) | n/a | Main text Eqs 6, 7 |
| Inhibition form E_max = E_max(0) * AR50 / (AR50 + AR_on) | n/a | Main text Eq 4 and Results paragraph on maximum inhibition |
| k_SR = ((kgrowth - kdeath) / bmax) * (S + R) * [S + R > BP] | n/a | Methods paragraph on transfer rate, with derivation of beta from Bmax |

## Units (dimensional analysis)

| Symbol | Meaning | Units |
|--------|---------|-------|
| `bact_growing, bact_resting` | susceptible / resting bacterial counts | CFU/mL |
| `ar_off, ar_on` | adaptive-resistance binding-model states | unitless (mass-balance 1) |
| `cgent` | gentamicin concentration in flask | mg/L |
| `kgrowth, kdeath, kel = exp(lkel), drug, emax, emax0, ksr` | rate constants | 1/h |
| `kon` | binding on-rate | L/(mg*h) |
| `koff` | binding off-rate | 1/h |
| `ec50` | half-effect concentration | mg/L |
| `ar50` | half-effect AR_on level | unitless |
| `bp, bmax, s0` | bacterial count scales | CFU/mL |
| `beta` | inverse-carrying-capacity constant | 1/(CFU/mL * h) |

Every bacterial ODE term has the form (1/h) * (CFU/mL) = (CFU/mL)/h; the
AR-binding ODEs are (L/(mg*h)) * (mg/L) * (unitless) = (unitless)/h; the
`cgent` ODE is (1/h) * (mg/L) = (mg/L)/h. All consistent.

```{r load-model}
mod <- rxode2::rxode(readModelDb("Mohamed_2012_gentamicin"))
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 / fixed parameter values as loaded from the model file (Mohamed 2012 Table 1)."
)
```

# Mechanistic checks

## 1. Maximum adaptive-resistance inhibition

The paper reports that "the maximum E_max inhibition by the adaptive resistance
mechanism was 90% as the maximum amount possible in AR_on is 1 [1/(1 + 0.113)]"
(Results paragraph below Table 1). Mass balance keeps `ar_on + ar_off = 1`, so
the saturating AR_on is 1 and the inhibition fraction is 1/(1 + AR50). With
AR50 = 0.113 this is 0.898 -- within rounding of the published 90%.

```{r ar-inhibition}
ar50 <- 0.113
data.frame(
  AR_on = c(0, 0.05, ar50, 0.5, 1),
  inhibition_fraction = c(0, 0.05, ar50, 0.5, 1) / (ar50 + c(0, 0.05, ar50, 0.5, 1)),
  E_max_remaining_fraction = ar50 / (ar50 + c(0, 0.05, ar50, 0.5, 1))
) |>
  knitr::kable(digits = 3,
               caption = "Adaptive-resistance attenuation of E_max. At AR_on = 1 the inhibition fraction is 0.898 (90%), matching Mohamed 2012.")
```

## 2. Adaptive-resistance half-lives at constant gentamicin

For first-order saturation of `ar_on` toward its steady state, the relaxation
rate is `kon * C + koff` and the half-time is `ln(2) / (kon * C + koff)`. The
paper reports concentration-dependent half-lives of 16, 4, and 1 h at constant
gentamicin concentrations of 1, 4, and 16 mg/L (Mohamed 2012 Discussion). The
binding-model arithmetic reproduces the 4-mg/L and 16-mg/L values closely; at
1 mg/L the simple two-state half-life is shorter (~12 h) than the published 16 h
estimate, reflecting that at low C the off-rate dominates and the published
"half-life" was likely read off the simulated trajectory rather than computed
analytically.

```{r ar-halflife}
kon  <- 0.0426
koff <- 0.0139
data.frame(
  C_mgL = c(1, 4, 16),
  rate_per_h = kon * c(1, 4, 16) + koff,
  half_life_h_model = log(2) / (kon * c(1, 4, 16) + koff),
  half_life_h_paper = c(16, 4, 1)
) |>
  knitr::kable(digits = 2,
               caption = "Adaptive-resistance half-life ln(2)/(kon*C + koff) versus the values cited in the Discussion of Mohamed 2012.")
```

# Growth-control replication (Fig 2A, antibiotic-free curve)

With no gentamicin, the population grows exponentially from the inoculum
(`s0 = 4.83e5`) and approaches the stationary `bmax = 8.26e8` CFU/mL plateau
once the total exceeds `bp = 2.09e6`. The paper reports a plateau "of
approximately 1e9 CFU/mL" in growth-control experiments -- consistent with
the model's `bmax` of 8.26e8.

```{r growth-control}
ev_gc <- rxode2::et(seq(0, 24, by = 0.25))
gc <- rxode2::rxSolve(mod, ev_gc, returnType = "data.frame", maxsteps = 1e5)

cat(sprintf("Inoculum  ln(CFU/mL)   = %.3f  (S(0) = %.2e)\n",
            gc$Cc[1], 4.83e5))
cat(sprintf("Plateau   ln(CFU/mL)   = %.3f  (B_max = 8.26e8, ln = %.3f)\n",
            tail(gc$Cc, 1), log(8.26e8)))

ggplot(gc, aes(time, Cc / log(10))) +
  geom_line(linewidth = 1) +
  geom_hline(yintercept = log10(8.26e8), linetype = 2, colour = "grey50") +
  labs(x = "Time (h)", y = expression(log[10]~CFU/mL),
       title = "Antibiotic-free growth control (replicates Fig 2A grey curve)",
       caption = "Dashed line: stationary B_max = 8.26e8 CFU/mL.")
```

# Static-experiment replication (Fig 2A)

Figure 2A shows seven constant gentamicin concentrations from 0.125 to 16 mg/L
plus a growth control. The paper text reports that "regrowth occurred for
concentrations of 0.5, 1, and 2 mg/L while the bacterial count remained below
the LOD during the 24 h of incubation following exposures to concentrations
of 4 and 16 mg/L" (Results, Time-kill curve experiments). Static experiments
correspond to `kel = 0` (gentamicin does not decline).

```{r fig2a-static}
loc <- 10            # limit of detection (CFU/mL)
conc <- c(0.5, 1, 2, 4, 16)

run_static <- function(C) {
  ev <- rxode2::et(amt = C, cmt = "cgent", time = 0) |>
    rxode2::et(seq(0, 24, by = 0.5))
  rxode2::rxSolve(mod, ev,
                  params = c(lkel = log(1e-12)),  # static: no decline
                  returnType = "data.frame", maxsteps = 1e5) |>
    mutate(C_mgL = C)
}

sim_static <- bind_rows(lapply(conc, run_static))

sim_static |>
  mutate(C_label = factor(sprintf("%g mg/L", C_mgL),
                          levels = sprintf("%g mg/L", conc))) |>
  mutate(log10_CFU = pmax(log10(pmax(bact_growing + bact_resting, 1)),
                          log10(loc))) |>
  ggplot(aes(time, log10_CFU, colour = C_label)) +
  geom_line(linewidth = 0.9) +
  geom_hline(yintercept = log10(loc), linetype = 3, colour = "grey50") +
  labs(x = "Time (h)", y = expression(log[10]~CFU/mL),
       colour = "Gentamicin", title = "Static time-kill curves (replicates Fig 2A)",
       caption = "Dotted line: limit of detection (10 CFU/mL).")
```

```{r fig2a-endpoints}
sim_static |>
  group_by(C_mgL) |>
  summarise(
    nadir_log10  = round(log10(min(pmax(bact_growing + bact_resting, 1))), 2),
    at_24h_log10 = round(log10(pmax(tail(bact_growing + bact_resting, 1), 1)), 2),
    .groups = "drop"
  ) |>
  knitr::kable(caption = "Nadir and 24-h total bacterial count by static gentamicin concentration. The model reproduces regrowth at 0.5-2 mg/L and below-LOD bacterial counts at 4 and 16 mg/L (Mohamed 2012 Fig 2A, Results paragraph 1).")
```

# Dynamic-experiment replication (Fig 2B single doses)

The dynamic experiments declined gentamicin biphasically: rate constant
0.33 /h for the first 4 h, then 0.037 /h thereafter. We simulate this by
splitting the run into two segments around the 4-h flow-rate switch.

```{r fig2b-dynamic}
peak <- c(2.0, 3.9, 7.8, 16)

run_dynamic <- function(Cmax) {
  ev1 <- rxode2::et(amt = Cmax, cmt = "cgent", time = 0) |>
    rxode2::et(seq(0, 4, by = 0.25))
  s1 <- rxode2::rxSolve(mod, ev1, params = c(lkel = log(0.33)),
                        returnType = "data.frame", maxsteps = 1e5)
  init2 <- tail(s1, 1)
  ev2 <- rxode2::et(seq(4, 24, by = 0.5))
  s2 <- rxode2::rxSolve(mod, ev2,
                        params = c(lkel = log(0.037)),
                        inits = c(bact_growing = init2$bact_growing,
                                  bact_resting = init2$bact_resting,
                                  ar_off       = init2$ar_off,
                                  ar_on        = init2$ar_on,
                                  cgent        = init2$cgent),
                        returnType = "data.frame", maxsteps = 1e5)
  bind_rows(s1, s2) |> mutate(Cmax_mgL = Cmax)
}

sim_dyn <- bind_rows(lapply(peak, run_dynamic))

p_left <- sim_dyn |>
  mutate(Cmax_label = factor(sprintf("Cmax %g mg/L", Cmax_mgL),
                             levels = sprintf("Cmax %g mg/L", peak))) |>
  ggplot(aes(time, cgent, colour = Cmax_label)) +
  geom_line(linewidth = 0.9) +
  scale_y_log10() +
  labs(x = "Time (h)", y = "Gentamicin (mg/L, log scale)",
       colour = NULL, title = "Simulated in-vitro PK (Fig 2B input)") +
  theme(legend.position = "bottom")

p_right <- sim_dyn |>
  mutate(Cmax_label = factor(sprintf("Cmax %g mg/L", Cmax_mgL),
                             levels = sprintf("Cmax %g mg/L", peak))) |>
  mutate(log10_CFU = pmax(log10(pmax(bact_growing + bact_resting, 1)),
                          log10(10))) |>
  ggplot(aes(time, log10_CFU, colour = Cmax_label)) +
  geom_line(linewidth = 0.9) +
  geom_hline(yintercept = log10(10), linetype = 3, colour = "grey50") +
  labs(x = "Time (h)", y = expression(log[10]~CFU/mL),
       colour = NULL, title = "Replicates Fig 2B (dynamic single dose)") +
  theme(legend.position = "bottom")

print(p_left)
print(p_right)
```

The model reproduces the paper's qualitative Fig 2B finding: regrowth
occurred for all dynamic single doses, but the higher peaks produced
deeper initial kills before regrowth.

```{r fig2b-endpoints}
sim_dyn |>
  group_by(Cmax_mgL) |>
  summarise(
    nadir_log10  = round(log10(min(pmax(bact_growing + bact_resting, 1))), 2),
    at_24h_log10 = round(log10(pmax(tail(bact_growing + bact_resting, 1), 1)), 2),
    ar_on_max    = round(max(ar_on), 3),
    .groups = "drop"
  ) |>
  knitr::kable(caption = "Nadir and 24-h total bacterial count plus peak AR_on after a single dynamic gentamicin dose. The 16 mg/L peak drives AR_on close to its saturation of 1, so a subsequent dose would face strongly reduced E_max (Mohamed 2012 Fig 2C resistance development).")
```

# Adaptive-resistance trajectory under a saturating dose

A single 16 mg/L peak (the highest dynamic experiment in Fig 2B) drives
`ar_on` past 0.5 within ~1 h, demonstrating that gentamicin rapidly induces
adaptive resistance even as it kills bacteria. After the gentamicin declines
(via the two-phase in-vitro PK), `ar_on` falls back through `koff` with a
characteristic half-life of ln(2)/0.0139 = 50 h (paper's fixed value).

```{r ar-trajectory}
ar_run <- sim_dyn |> filter(Cmax_mgL == 16)
ggplot(ar_run, aes(time)) +
  geom_line(aes(y = ar_on, colour = "ar_on"), linewidth = 1) +
  geom_line(aes(y = cgent / 16, colour = "cgent / 16 (mg/L)"), linewidth = 0.6, linetype = 2) +
  labs(x = "Time (h)", y = "AR_on (solid) or scaled cgent (dashed)",
       colour = NULL,
       title = "Adaptive-resistance development under a 16 mg/L peak (Fig 2B context)")
```

# Assumptions and deviations

* **Model class / species.** This is an in-vitro semi-mechanistic PKPD model, not
  a popPK model; `population$species` records the *E. coli* ATCC 25922 isolate.
  No PKNCA validation is performed (there is no drug NCA to compute); the
  mechanistic checks above replace it.
* **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 gentamicin, so the model file and
  this vignette use `Mohamed_2012_gentamicin`.
* **Final model.** Parameters are the final population-mean estimates from
  Table 1; alternative AR formulations explored in the paper (time-and-concn
  EC50 of Tam et al. 2005, turnover EC50, third bacterial compartment) all gave
  higher OFV and are not packaged.
* **Fixed parameters.** `kdeath = 0.179 /h` is carried over from Nielsen 2007
  (ref 36 in the paper). `koff = 0.0139 /h` was fixed at the lowest value that
  did not worsen OFV (the data only weakly identify the return-to-susceptibility
  rate). Both are wrapped in `fixed()` in `ini()`. `lkel = log(0.037)` is the
  default (dynamic phase-2, the neonatal terminal kel); pass
  `params = c(lkel = log(1e-12))` to replicate static experiments (no decline)
  or override to `log(0.33)` then `log(0.037)` to replicate the dynamic
  in-vitro kinetic system (see the Fig 2B chunk).
* **In-vitro flask elimination.** The gentamicin concentration is a state
  variable; the user adds gentamicin via dosing events. For dynamic experiments
  the paper used two flow-rate phases (0.33 /h for the first 4 h, then 0.037 /h
  thereafter) to mimic preterm-neonate kinetics. We replicate this by splitting
  the simulation into two segments rather than coding a piecewise `kel` into
  the ODE, so the model file stays general for any user-supplied PK driver.
* **Initial inoculum.** `s0 = 4.83e5 CFU/mL` is the average across all
  experimental start inocula (Methods) and the value the paper used for its
  dosing-schedule predictions. The high-inoculum subset (12 h pre-growth to
  ~1e9 CFU/mL) used to study the inoculum effect is not packaged as a separate
  scenario; users can simulate it by changing `s0` or by initializing
  `bact_resting` non-zero.
* **Three residual-error magnitudes were estimated.** Mohamed 2012 fit
  `RE_static = 1.69`, `RE_dynamic = 2.80`, and a replicate-specific
  `RRE = 0.618` (all on the natural-log scale; Methods + Table 1). The model
  file ships `addSd = 1.69` (RE_static) as the single residual error because
  the conditional `RE_static | RE_dynamic` switch would require a new
  experiment-type covariate. For VPC-style simulation against dynamic
  experiments, scale `addSd` to 2.80 at rxSolve time.
* **No inter-experiment variability.** The paper did not estimate
  inter-experiment variability ("interexperimental variability was not
  estimated"); accordingly the model file has no `eta*` IIV terms.
* **k_SR breakpoint discontinuity.** Mohamed defines `k_SR` as zero below the
  estimated bacterial-count breakpoint `BP` and a linear function above it.
  This is encoded with the indicator `(total_bact > bp)`, which is a hard
  discontinuity in the right-hand side. Pass `maxsteps = 1e5` to `rxSolve`
  for robust integration through the threshold; the bundled simulations all
  do.
* **High-inoculum-effect overprediction.** For the auxiliary 12-h-grown
  high-inoculum experiments (Fig 3C in the paper), the model overpredicts
  killing at 2 and 4 mg/L. This is a documented model limitation -- the model
  does not represent the inoculum effect on susceptibility; the high-inoculum
  subset is not validated in this vignette.
* **Upstream popPK is out of scope.** The paper's neonatal dosing-schedule
  predictions (Figs 4-6) linked this PD model to a separate 3-compartment
  popPK developed in Nielsen 2011 (ref 35). That popPK is not packaged here.
  Users wishing to reproduce the neonatal predictions should drive `cgent`
  with the upstream popPK output.
* **Convention deviations** (`checkModelConventions()` warnings, no errors).
  All are expected for an in-vitro mechanism-based PKPD model: (a) the
  bacterial-state (`bact_growing`, `bact_resting`), adaptive-resistance
  (`ar_off`, `ar_on`), and gentamicin-concentration (`cgent`) compartments
  are paper-mechanistic, declared via `paper_specific_compartment_pattern`;
  (b) the single observation `Cc` carries a non-PK output (natural log of
  total viable count, not a drug concentration); (c) the dosing/concentration
  units are both `mg/L` because the antibiotic input is a concentration in
  the in-vitro system.
