---
title: "Colistin meropenem (Mohamed 2016)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Colistin meropenem (Mohamed 2016)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: Mohamed A. F., Kristoffersson A. N., Karvanen M., Nielsen E. I., Cars O., Friberg L. E. (2016). Dynamic interaction of colistin and meropenem on a WT and a resistant strain of *Pseudomonas aeruginosa* as quantified in a PK/PD model. *J Antimicrob Chemother* 71(5):1279-1290. doi:10.1093/jac/dkv488.
* DDMORE Foundation Model Repository: [DDMODEL00000173](https://repository.ddmore.eu/model/DDMODEL00000173).
* Article: <https://doi.org/10.1093/jac/dkv488>.

The model is a 10-compartment in vitro time-kill PK/PD model that quantifies
the joint pharmacodynamic effect of colistin and meropenem on a wild-type
P. aeruginosa strain (ATCC 27853) and a meropenem-resistant clinical isolate
(ARU552). The bacterial system has proliferating (S) and resting (R)
subpopulations linked through a Bmax-driven feedback term, plus a colistin
adaptive-resistance compartment pair (`col_ce` / `col_bindoff`) and a
meropenem-resistant mutant subpopulation (`S_mut` / `R_mut`). Drug
concentrations decline by first-order rate constants `KE1` (meropenem,
constant 0.02 1/h) and `KE2` (colistin, a 50+ entry per-experiment lookup
keyed on the initial colistin concentration `Ccol`).

# Population

The dataset is in vitro bacterial culture, not human subjects. Two strains
are studied:

* **BACT = 2** — *P. aeruginosa* ATCC 27853 (wild-type, drug-susceptible).
* **BACT = 1** — *P. aeruginosa* ARU552 (meropenem-resistant clinical
  isolate from a Swedish ICU patient).

Experimental conditions (`TYPE`):

* `TYPE = 1` meropenem alone, `TYPE = 2` colistin alone, `TYPE = 3`
  combination, `TYPE = 4` / `TYPE = 5` no-drug control.

Static and dynamic time-kill experiments were combined; `ID2 > 18` flags
dynamic experiments where bacterial growth rate is scaled by
`1 + THETA(35)`. Initial inocula range from 2.26e5 to 3.5e5 CFU/mL across
experiments. Replicate-specific residual error (`REPL`) and below-LOQ
censoring (`BLOQ`, M3 method, LLOQ = 10 CFU/mL) are handled by NONMEM in
the source; the nlmixr2 extraction collapses replicates into a single
combined SD and relies on the standard `cens` column for BLOQ.

The full population description is available at
`readModelDb("Mohamed_2016_colistin_meropenem")$population`.

# Source trace

Every parameter line in
`inst/modeldb/ddmore/Mohamed_2016_colistin_meropenem.R` carries a trailing
comment pointing to the THETA index in
`Output_real_ColistinMeropenem_interaction.lst` and the source line in
`Executable_ColistinMeropenem_Interaction.mod`. The block below reproduces
the cross-walk for review.

| nlmixr2 parameter           | NONMEM THETA | Final value (.lst, .mod $THETA initial) | Source location |
|-----------------------------|--------------|-----------------------------------------|-----------------|
| `kga`                       | THETA(1)     | 1.08 1/h                                | .mod L485, .lst FINAL PARAMETER ESTIMATE block (line 855) |
| `kgp`                       | THETA(2)     | 0.814 1/h                               | .mod L486, .lst L855 |
| `kk`                        | THETA(3) FIX | 0.179 1/h                               | .mod L487, .lst L855 |
| `emax_col_atcc`             | THETA(4) FIX | 282                                     | .mod L488, .lst L855 |
| `slop1_col_aru`             | THETA(5)     | 17.2                                    | .mod L489, .lst L855 |
| `emax_mero_aru`             | THETA(7)     | 1.74 1/h                                | .mod L491, .lst L855 |
| `ec51_col_atcc`             | THETA(8)     | 1.84 mg/L                               | .mod L492, .lst L855 |
| `slopam_mero_atcc`          | THETA(9)     | 2.16                                    | .mod L493, .lst L855 |
| `ec53_mero_aru`             | THETA(10)    | 17.7 mg/L                               | .mod L494, .lst L855 |
| `gamc_col`                  | THETA(11) FIX| 1                                       | .mod L495, .lst L855 |
| `gama_mero_atcc`            | THETA(12)    | 0.376                                   | .mod L496, .lst L855 |
| `bmax_log10`                | THETA(13)    | 9.11 (BMAX = 10^9.11 ≈ 1.29e9 CFU/mL)   | .mod L497, .lst L856 |
| `kon51_col_atcc`            | THETA(14)    | 2.42 mg/L                               | .mod L498, .lst L856 |
| `kon52_col_aru`             | THETA(15)    | 18.6 mg/L                               | .mod L499, .lst L856 |
| `koffc_col`                 | THETA(16)    | 0.103 1/h                               | .mod L500, .lst L856 |
| `koffm_mero` (FIX 0)        | THETA(17) FIX| 0                                       | .mod L501, .lst L856 |
| `konca_col_atcc`            | THETA(18)    | 0.779 1/h                               | .mod L502, .lst L856 |
| `koncp_col_aru`             | THETA(19)    | 2.68 1/h                                | .mod L503, .lst L856 |
| `konma_mero_atcc` (FIX)     | THETA(20) FIX| 0.113                                   | .mod L504, .lst L856 |
| `konmp_mero_aru` (FIX)      | THETA(21) FIX| 0.0821                                  | .mod L505, .lst L856 |
| `gam3c_col`                 | THETA(22)    | 1.08                                    | .mod L506, .lst L856 |
| `gam1_col`                  | THETA(23)    | 0.0245                                  | .mod L507, .lst L856 |
| `gamp_mero_aru`             | THETA(26)    | 2.79                                    | .mod L510, .lst L857 |
| `inter_atcc`                | THETA(27)    | -2.58                                   | .mod L511, .lst L857 |
| `inter_aru`                 | THETA(28)    | 1.30                                    | .mod L512, .lst L857 |
| `neg_log10_mutf_atcc`       | THETA(29)    | 2.83 (MUTA = 10^-2.83 ≈ 1.48e-3)        | .mod L513, .lst L857 |
| `concshift_mut`             | THETA(31)    | 22.5                                    | .mod L515, .lst L857 |
| `kg_decrease_atcc` (FIX 0)  | THETA(33) FIX| 0                                       | .mod L517, .lst L857 |
| `kg_decrease_aru`           | THETA(34)    | 0.536                                   | .mod L518, .lst L857 |
| `kg_dyn_increase`           | THETA(35)    | 0.332                                   | .mod L519, .lst L857 |
| `addSd_log_cfu`             | SIGMA(1+2)   | sqrt(0.422 + 0.0364) * log(10) ≈ 1.559  | .mod L521-526, .lst SIGMA block lines 877-889 |

NONMEM placeholders THETA(6) (FIX 0), THETA(24) (FIX 1), THETA(25) (FIX 1),
THETA(30) (FIX 0), and THETA(32) (FIX 0) are not referenced anywhere in
`$DES` / `$ERROR`; they are deliberately omitted from `ini()` because
nlmixr2 requires every `ini()` entry to be used inside `model()`. The
omission is documented inline in the model file's `ini()` block.

The ODE structure is reproduced verbatim from `$DES` (lines 280-437 of the
.mod). The 56-line KE2 lookup cascade in `$PK` (lines 199-276) is
preserved verbatim — KE2 is set by exact-equality tests on the
experimental-well initial colistin concentration `Ccol`. Any `Ccol` value
not in the cascade leaves `KE2 = 0`.

# Validation strategy

The bundle ships an in vitro experimental dataset that is not a
re-extractable population PK study, and the linked publication is not on
local disk (searched
`/home/bill/github/mab_human_consensus/literature/`). The validation
strategy is therefore **F.2 self-consistency** per
`extract-literature-model` `references/verification-checklist.md`:
re-simulate the bundle's `Simulated_ColistinMeropenem_Interaction.csv`
through `rxode2::rxSolve()` with typical-value parameters (no IIV, no
residual error) and compare the trajectory against the simulated
observations in the bundle. PKNCA does not apply (this is bacterial-count
PD on log scale, not concentration-time PK). Population NCA comparison
against a published table is not possible because the publication is not
on disk.

# Re-simulation of the bundle's simulated dataset

```{r load-simulated-data}
sim_csv <- system.file("..", "vignettes", "articles", "data",
                       "Mohamed_2016_simulated.csv",
                       package = "nlmixr2lib", mustWork = FALSE)
if (!file.exists(sim_csv)) {
  sim_csv <- file.path("data", "Mohamed_2016_simulated.csv")
}
raw <- read.csv(sim_csv, stringsAsFactors = FALSE)
cat("Rows:", nrow(raw), "  unique IDs:", length(unique(raw$ID)), "\n")
cat("TYPE distribution:\n")
print(table(raw$TYPE, dnn = "TYPE"))
cat("BACT distribution:\n")
print(table(raw$BACT, dnn = "BACT"))
```

The CSV has 2,636 observation rows over 139 unique experiments. To keep
the vignette inside the 5-minute render gate while still covering every
experimental condition, we subset to one representative experiment per
(BACT × TYPE) cell with a non-trivial Ccol or Cmer.

```{r select-cohort}
# One representative experiment per (BACT, TYPE) cell. We prefer
# experiments whose Ccol matches a value in the .mod KE2 cascade so the
# colistin-elimination rate is non-zero where the source intends it.
representative_ids <- raw |>
  dplyr::filter(EVID == 1) |>
  dplyr::group_by(BACT, TYPE) |>
  dplyr::slice(1) |>
  dplyr::ungroup() |>
  dplyr::pull(ID)
cat("Representative IDs (one per (BACT, TYPE) cell):\n")
print(representative_ids)

events_raw <- raw |>
  dplyr::filter(ID %in% representative_ids) |>
  dplyr::select(ID, TIME, EVID, CMT, AMT, DV, Ccol, Cmer, TYPE, BACT,
                DILcol, DILmer, ID2, REPL, BLOQ) |>
  # The model expects the canonical CONMED_<drug>_CC covariate names
  # (`CONMED_COL_CC` for colistin, `CONMED_MER_CC` for meropenem).
  # Rename the source-dataset shorthand columns here so the events
  # table satisfies the model's covariate declarations.
  dplyr::rename(CONMED_COL_CC = Ccol, CONMED_MER_CC = Cmer)
cat("Filtered rows:", nrow(events_raw), "\n")
```

```{r build-events}
# Add a second EVID=1 dose to the S_mut compartment (CMT=9) with
# AMT = inoculum * MUT, mirroring the .mod $PK A_0(9) = NMUT
# initialization (NMUT = AMT * MUT, where MUT = 10^-neg_log10_mutf_atcc
# when MERO == 1, i.e. TYPE in {1, 3}; otherwise MUT = 0).
mut_atcc <- 10^(-2.83)        # MUTA per .mod L147
# .mod L148 sets MUTP = MUTA, so the strain split does not matter for the
# initial mutant fraction.
mut_for_row <- function(TYPE) ifelse(TYPE %in% c(1, 3), mut_atcc, 0)

dose_S <- events_raw |>
  dplyr::filter(EVID == 1) |>
  dplyr::mutate(CMT = 1L)           # S compartment

dose_Smut <- events_raw |>
  dplyr::filter(EVID == 1) |>
  dplyr::mutate(
    CMT = 9L,                                # S_mut compartment
    AMT = AMT * mut_for_row(TYPE)
  )

obs <- events_raw |>
  dplyr::filter(EVID == 0) |>
  dplyr::mutate(CMT = NA_integer_)

events <- dplyr::bind_rows(dose_S, dose_Smut, obs) |>
  dplyr::arrange(ID, TIME, dplyr::desc(EVID))

cat("Event rows after augmentation:", nrow(events),
    " (dose-S:", sum(events$EVID == 1 & events$CMT == 1),
    ", dose-Smut:", sum(events$EVID == 1 & events$CMT == 9),
    ", obs:", sum(events$EVID == 0), ")\n")
```

```{r typical-value-simulation}
mod <- readModelDb("Mohamed_2016_colistin_meropenem")
mod_typical <- rxode2::zeroRe(mod)
sim <- rxode2::rxSolve(mod_typical, events = events,
                       keep = c("DV", "BLOQ", "REPL", "TYPE", "BACT"),
                       returnType = "data.frame")
sim_df <- sim |>
  dplyr::filter(!is.na(DV) & DV > 0) |>      # keep observation rows
  dplyr::rename(time = time)
cat("Simulated rows:", nrow(sim_df), "\n")
```

# Self-consistency check

The DV column in the simulated CSV is a NONMEM-generated stochastic
observation (IPRED + per-replicate residual). Our typical-value
re-simulation produces `log_cfu` = ln(total bacteria) which corresponds
directly to NONMEM's IPRED (`$ERROR` line `IPRED=LOG(ATOT)`). For each
observation we compare `log_cfu` (rxode2 typical-value) against `DV`
(NONMEM-simulated stochastic observation).

```{r consistency-summary}
diff_df <- sim_df |>
  dplyr::mutate(
    cohort = paste0("BACT=", BACT, " / TYPE=", TYPE),
    residual_log = DV - log_cfu,
    abs_residual = abs(residual_log)
  )

summary_df <- diff_df |>
  dplyr::group_by(cohort) |>
  dplyr::summarise(
    n              = dplyr::n(),
    median_dv      = stats::median(DV),
    median_log_cfu = stats::median(log_cfu),
    rmse           = sqrt(mean(residual_log^2)),
    median_abs_res = stats::median(abs_residual),
    .groups = "drop"
  )
knitr::kable(summary_df, digits = 3,
             caption = paste("Per-cohort self-consistency summary.",
                             "RMSE / median |DV - log_cfu| are on natural-log scale."))
```

The expected per-observation residual under perfect translation is the
NONMEM stochastic noise term, with theoretical SD ≈
`sqrt(0.422 + 0.0364) * log(10)` ≈ `1.56` on the natural-log scale. RMSE
values comparable to this magnitude indicate the trajectory is consistent
with the source; values much larger flag a translation defect.

```{r consistency-figure, fig.cap = "Typical-value rxode2 trajectory (line) vs NONMEM-simulated DV (points), one panel per (BACT, TYPE) cohort."}
ggplot(diff_df, aes(time, log_cfu)) +
  geom_line(color = "steelblue") +
  geom_point(aes(y = DV), alpha = 0.5, size = 1) +
  facet_wrap(~ cohort, scales = "free_y") +
  labs(x = "Time (h)", y = "ln(total bacterial count, CFU/mL)",
       title = "Self-consistency: rxode2 typical value vs NONMEM-simulated DV",
       caption = paste("Points: DV column (NONMEM stochastic);",
                       "line: rxode2::rxSolve typical-value `log_cfu`."))
```

# Assumptions and deviations

The verbatim translation deliberately preserves several non-standard
features of the source model. Each is flagged here for downstream review.

* **Non-canonical observation variable.** `log_cfu` (ln of total
  bacterial count) replaces the canonical PK observation `Cc`, mirroring
  the in vitro PD precedent set by
  `inst/modeldb/therapeuticArea/oncology/oncology_sdm_lobo_2002.R`
  (which uses `tumorVol`). Flagged as a `checkModelConventions()`
  warning.
* **Non-canonical compartment names.** `S` / `R` / `S_mut` / `R_mut`
  (bacterial subpopulations); `mero` / `mero_ce` / `mero_bindoff` (the
  meropenem concentration plus two declared-but-unused adaptive-resistance
  states); `col` / `col_ce` / `col_bindoff` (colistin concentration plus
  active adaptive-resistance pair). None of these match the canonical
  `depot` / `central` / `peripheral*` / `effect` / `target` / `complex`
  list. The names follow the .mod $MODEL declarations directly. Flagged
  as 10 `checkModelConventions()` compartment warnings.
* **Dataset-tied covariates.** `BACT`, `TYPE`, `Ccol`, `Cmer`, `DILcol`,
  `DILmer`, and `ID2` are not in
  `inst/references/covariate-columns.md`. They are intrinsic to the in
  vitro experimental design (strain index, condition index, initial
  drug-well concentrations, bath dilution rates, within-subject experiment
  index for dynamic-system flagging) rather than reusable
  human-population covariates, so they are not promoted to the canonical
  register. Flagged as 7 `checkModelConventions()` covariate warnings.
* **Dosing / concentration unit dimensions intentionally differ.** The
  inoculum dose is in CFU/mL (a count density) while drug concentration
  units inside `mero` / `col` are in mg/L. There is no shared dimension
  to reconcile because the two compartment classes describe different
  physical quantities. Flagged as 1 `checkModelConventions()` units
  warning.
* **`KE2` cascade preserved verbatim.** The 56-line per-experiment
  `KE2` lookup is reproduced exactly with floating-point equality tests
  on `Ccol`. Any `Ccol` value not in the cascade leaves `KE2 = 0`.
  Downstream users repurposing the model with arbitrary colistin
  concentrations will need to set their own `KE2` (e.g., by editing
  the `model()` block or supplying `KE2` as a covariate). The verbatim
  decision was made by the operator
  (`queue/.claude_task_runner/sidecar/006-mohamed_2016_colistin_meropenem/response-001.json`).
* **`ETA(1)` dropped.** The .mod declares `ETA(1)` shared between `KGA`
  and `KGP` with `OMEGA = 0 FIX`, contributing nothing to the
  predictions. The IIV term is dropped from the nlmixr2 translation; the
  growth rates are typical values only.
* **Per-replicate residual error collapsed.** The .mod's `$ERROR` block
  uses one common epsilon plus four replicate-specific epsilons (all
  `BLOCK(1) SAME`, all `FIXED` at the same variance). The verbatim
  translation collapses these into a single combined additive residual
  on natural-log scale (`addSd_log_cfu = sqrt(0.422 + 0.0364) * log(10)`).
* **No linked publication on disk.** The associated paper
  ([doi:10.1093/jac/dkv488](https://doi.org/10.1093/jac/dkv488)) is not
  available locally; PKNCA-style comparison against published tables is
  therefore not possible. Validation reduces to F.2 self-consistency
  against the bundle-shipped `Output_simulated_*.lst`.
* **Unused THETA placeholders.** THETA(6), THETA(24), THETA(25),
  THETA(30), and THETA(32) are FIX placeholders in the source that are
  declared but never referenced in any `$PK` / `$DES` / `$ERROR`
  expression. They are documented inline as comments in the model
  file's `ini()` block but not declared as `ini()` parameters because
  nlmixr2 requires every `ini()` entry to appear in `model()`.
* **Mutant initial condition sourced via vignette preprocessing.** The
  .mod's `$PK` block sets `A_0(9) = NMUT = AMT * MUT` per dose record.
  rxode2 cannot read `AMT` from the event table inside the model body,
  so the vignette injects a second `EVID = 1` dose to compartment 9
  (`S_mut`) with the appropriate `AMT` per row. This produces the same
  initial condition for the mutant subpopulation. Downstream users
  applying the model to new datasets must replicate this preprocessing
  step.
