---
title: "Colistin (Leuppi-Taegtmeyer 2019)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Colistin (Leuppi-Taegtmeyer 2019)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: Leuppi-Taegtmeyer A, Decosterd LA, Osthoff M, Mueller NJ, Buclin T,
  Corti N. Multicenter Population Pharmacokinetic Study of Colistimethate Sodium
  and Colistin Dosed as in Normal Renal Function in Patients on Continuous Renal
  Replacement Therapy. Antimicrob Agents Chemother. 2019;63(2):e01957-18.
* Article: <https://doi.org/10.1128/AAC.01957-18>
* DDMORE Foundation Model Repository entry: DDMODEL00000295
  (<https://repository.ddmore.eu/model/DDMODEL00000295>).
* Local DDMORE bundle used for this extraction:
  `dpastoor/ddmore_scraping/295/` (Executable: `Executable_CMS_colistin_PK_CRRT.mod`,
  Output: `Output_real_CMS_colistin_PK_CRRT.lst`, simulated dataset:
  `Simulated_Data_CMS_colistin_PK_CRRT.csv`).

The published article was NOT on disk during extraction; all parameter values
and equations come from the DDMORE bundle's `Output_real_*.lst` (final
estimates) and `Executable_*.mod` (structural equations). The
`Model_Accommodations.txt` field "There are no model differences with the
publication referenced" indicates that the bundle reproduces the published
final-model structure verbatim. The vignette's validation strategy is therefore
F.2 self-consistency (re-simulate against the bundle's `Simulated_*.csv`)
rather than a publication side-by-side comparison.

# Population

The DDMORE bundle's `.rdf` model description states: "CMS and colistin
pharmacokinetics were assessed prospectively in 10 critically ill patients
requiring CRRT. Extensive pharmacokinetic sampling was performed on treatment
day 1, 3 and 5 after administration of a loading dose of CMS, followed by a
maintenance dosage every eight hours." Detailed baseline demographics (age,
weight, sex, study sites, comorbidities) live in the linked publication and
were not available during this extraction.

The bundle's simulated dataset
(`Simulated_Data_CMS_colistin_PK_CRRT.csv`, 11 subjects on a 7-day regimen)
reports per-subject covariate values:

* Age: 53 to 79 years (median 68);
* Body weight: 60 to 92 kg (median 75); SEX recorded as 1 or 2 (encoding not
  documented in the bundle).
* CRRT operational settings:
  - Blood flow `QBL` 100-150 mL/min (most subjects 100 or 150);
  - Effluent flow `QEFF` 1700-3000 mL/h;
  - Hematocrit `HT` is missing (`.`) for every subject in the simulated CSV;
    the .mod treats `HT == 0` as the missing-value sentinel and substitutes
    the population-typical fraction of 0.25.
* Dose pattern: loading dose 9 MIU CMS (720 mg of CMS sodium, ~300 mg colistin
  base equivalents at 1 MIU = 80 mg CMS Na = ~33.3 mg colistin base) infused
  over ~26 minutes, followed by 3 MIU q8h maintenance (240 mg CMS sodium per
  dose) for the remainder of the observation window.

The same information is available programmatically via the model's
`population` metadata
(`readModelDb("Leuppi-Taegtmeyer_2019_colistin")$population`).

# Source trace

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

| Equation / parameter        | Value | Source location                                                           |
|-----------------------------|-------|---------------------------------------------------------------------------|
| `lcl` = log(CL1M)           | log(2.31)  | `Output_real_*.lst` FINAL THETA TH 1; `.mod` line 116 (`(0, 2.31)`)  |
| `lvc` = log(V1)             | log(12.1)  | `Output_real_*.lst` FINAL THETA TH 2; `.mod` line 117                |
| `sc_cms` = SIV3             | 1.05       | `Output_real_*.lst` FINAL THETA TH 3; `.mod` line 118                |
| `lcl_col` = log(CL2M)       | log(1.93)  | `Output_real_*.lst` FINAL THETA TH 4; `.mod` line 119                |
| `lvc_col` = log(V2)         | log(70.1)  | `Output_real_*.lst` FINAL THETA TH 5; `.mod` line 120                |
| `sc_col` = SIV4             | 0.454      | `Output_real_*.lst` FINAL THETA TH 6; `.mod` line 121                |
| `propSd` (Cc and Cc_col)    | 0.222      | `Output_real_*.lst` FINAL THETA TH 7; `.mod` line 122                |
| `addSd`  (Cc and Cc_col)    | 0.459      | `Output_real_*.lst` FINAL THETA TH 8; `.mod` line 123                |
| OMEGA(1,1) = var(eta lcl)        | 0.270  | `Output_real_*.lst` FINAL OMEGA(1,1)                                |
| OMEGA(2,2) = var(eta lvc)        | 0.133  | `Output_real_*.lst` FINAL OMEGA(2,2)                                |
| OMEGA(3,3) = var(eta lcl_col)    | 0.304  | `Output_real_*.lst` FINAL OMEGA(3,3)                                |
| OMEGA(4,4) = var(eta lvc_col)    | 0.247  | `Output_real_*.lst` FINAL OMEGA(4,4)                                |
| Bioavailability factor F1   | 1155.5/1749.8 | `.mod` line 76 (mass-conversion factor: mg CMS Na -> mg colistin base) |
| Filter priming volume       | 0.2 L      | `.mod` line 34 (`V3 = 0.2`) and line 52 (`V4 = 0.2`)                 |
| Cartridge priming volume    | 0.3 L      | `.mod` line 43 (`V5 = 0.3`) and line 61 (`V6 = 0.3`)                 |
| ODE: `d/dt(central)`        | n/a        | `.mod` `$DES` line 87 (DADT(1))                                      |
| ODE: `d/dt(central_col)`    | n/a        | `.mod` `$DES` line 88 (DADT(2))                                      |
| ODE: `d/dt(filter)`         | n/a        | `.mod` `$DES` line 89 (DADT(3))                                      |
| ODE: `d/dt(filter_col)`     | n/a        | `.mod` `$DES` line 90 (DADT(4))                                      |
| ODE: `d/dt(cart)`           | n/a        | `.mod` `$DES` line 91 (DADT(5))                                      |
| ODE: `d/dt(cart_col)`       | n/a        | `.mod` `$DES` line 92 (DADT(6))                                      |
| ODE: `d/dt(effl)`           | n/a        | `.mod` `$DES` line 93 (DADT(7))                                      |
| ODE: `d/dt(effl_col)`       | n/a        | `.mod` `$DES` line 94 (DADT(8))                                      |

NONMEM cumulative-AUC compartments DADT(9) / DADT(10) are dropped because PKNCA
in this vignette computes AUC analytically.

# Virtual cohort

The original observed data are not publicly available. The simulations below
use a virtual cohort that mimics the bundle's `Simulated_Data_*.csv` regimen:
a 9 MIU loading dose (720 mg CMS sodium) over 26 minutes followed by 3 MIU
maintenance (240 mg CMS sodium) infused over 60 minutes every 8 h to t = 168 h
(7 days). CRRT settings are set to the median of the bundle (QBL = 100 mL/min,
QEFF = 2100 mL/h) with hematocrit at the model's default 0.25 (the bundle's
simulated subjects all carry `HT = .`).

```{r cohort}
set.seed(20250506)

dose_mg_cms_na  <- 720    # loading dose, mg CMS sodium
maint_mg_cms_na <- 240    # maintenance dose, mg CMS sodium
load_inf_h      <- 26 / 60   # 26-min infusion for the loading dose
maint_inf_h     <- 60 / 60   # 60-min infusion for maintenance
sim_hours       <- 168      # 7 days

dose_times <- c(0, seq(8, sim_hours - 8, by = 8))
n_subjects <- 10L

build_subject_events <- function(id, qbl, qeff, hct) {
  # Dosing rows: loading + maintenance, all into central (CMS).
  doses <- data.frame(
    id   = id,
    time = dose_times,
    amt  = c(dose_mg_cms_na, rep(maint_mg_cms_na, length(dose_times) - 1L)),
    cmt  = "central",
    evid = 1L,
    rate = c(dose_mg_cms_na / load_inf_h,
             rep(maint_mg_cms_na / maint_inf_h, length(dose_times) - 1L))
  )
  # Observation grid: single observation cmt per time (rxode2 returns BOTH
  # Cc and Cc_col at every observation row, so duplicating the row per
  # output cmt would only produce duplicate (id, time) pairs that PKNCA
  # rejects).
  tgrid <- seq(0.001, sim_hours, by = 0.5)
  obs <- data.frame(
    id   = id,
    time = tgrid,
    cmt  = "Cc",
    amt  = NA_real_,
    rate = NA_real_,
    evid = 0L
  )
  ev <- dplyr::bind_rows(doses, obs)
  ev$QBL  <- qbl
  ev$QEFF <- qeff
  ev$HCT  <- hct
  ev[order(ev$id, ev$time, ev$evid), ]
}

cohort <- lapply(seq_len(n_subjects), function(i) {
  build_subject_events(id = i, qbl = 100, qeff = 2100, hct = 0.25)
}) |> dplyr::bind_rows()

stopifnot(!anyDuplicated(unique(cohort[, c("id", "time", "evid")])))
```

# Simulation

```{r simulate}
mod <- rxode2::rxode(readModelDb("Leuppi-Taegtmeyer_2019_colistin"))
sim <- rxode2::rxSolve(mod, events = cohort,
                       keep = c("QBL", "QEFF", "HCT")) |>
  as.data.frame()
```

For the figures we also produce a typical-value (no-IIV / no-residual-error)
trajectory using `rxode2::zeroRe()`:

```{r simulate-typical}
mod_typical <- mod |> rxode2::zeroRe()
sim_typical <- rxode2::rxSolve(mod_typical, events = cohort[cohort$id == 1L, ]) |>
  as.data.frame()
```

# Replicate published trajectories

The publication is not on disk; the figures below are drawn against the
DDMORE bundle's behaviour rather than a paper figure number. The plotted
ranges reproduce the qualitative shape of the simulated dataset shipped with
the bundle (approximately 30-50 mg/L peak CMS during the loading infusion;
approximately 2-5 mg/L colistin at steady state).

```{r figure-cms-trajectory}
sim |>
  dplyr::filter(time > 0, time <= sim_hours) |>
  dplyr::group_by(time) |>
  dplyr::summarise(
    Q05 = quantile(Cc, 0.05, na.rm = TRUE),
    Q50 = quantile(Cc, 0.50, na.rm = TRUE),
    Q95 = quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  ggplot(aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25) +
  geom_line(linewidth = 0.8) +
  labs(x = "Time (h)",
       y = "CMS plasma concentration (mg/L colistin equivalents)",
       title = "Stochastic CMS plasma trajectory (10-subject virtual cohort)",
       caption = "9 MIU loading + 3 MIU q8h IV. CRRT QBL = 100 mL/min, QEFF = 2100 mL/h, HCT = 0.25.") +
  scale_y_log10()
```

```{r figure-col-trajectory}
sim |>
  dplyr::filter(time > 0, time <= sim_hours) |>
  dplyr::group_by(time) |>
  dplyr::summarise(
    Q05 = quantile(Cc_col, 0.05, na.rm = TRUE),
    Q50 = quantile(Cc_col, 0.50, na.rm = TRUE),
    Q95 = quantile(Cc_col, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  ggplot(aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25) +
  geom_line(linewidth = 0.8) +
  labs(x = "Time (h)",
       y = "Colistin plasma concentration (mg/L)",
       title = "Stochastic colistin plasma trajectory (10-subject virtual cohort)",
       caption = "9 MIU loading + 3 MIU q8h IV. CRRT QBL = 100 mL/min, QEFF = 2100 mL/h, HCT = 0.25.")
```

# F.2 self-consistency: re-simulate the DDMORE-bundle scenario

Per `references/verification-checklist.md` Section F.2, the substitute for the
publication-comparison check is to reproduce the typical-value trajectory at
the bundle's reported regimen and confirm the shape of the trajectory.

```{r selfconsistency}
typical_summary <- sim_typical |>
  dplyr::filter(time > 0) |>
  dplyr::summarise(
    cms_max         = max(Cc, na.rm = TRUE),
    cms_max_time_h  = time[which.max(Cc)],
    col_steady_24h  = Cc_col[which.min(abs(time - 24))],
    col_steady_72h  = Cc_col[which.min(abs(time - 72))],
    col_steady_168h = Cc_col[which.min(abs(time - 168))],
    col_max         = max(Cc_col, na.rm = TRUE)
  )
knitr::kable(typical_summary, digits = 2,
             caption = "Typical-value simulation summary at the DDMORE-bundle regimen.")
```

The CMS peak occurs at the end of the loading infusion (t ~= 0.43 h, ~35 mg/L
colistin equivalents); the colistin trajectory rises gradually through the
first 24 h and approaches a plateau of approximately 4 mg/L by 72 h, in the
same magnitude as the simulated dataset's mid-day-3 colistin DV values
(approximately 2 to 5 mg/L) shipped with the bundle.

# PKNCA validation

PKNCA is run on the colistin trajectory (the active species and the clinical
PK target). Because the source publication is not on disk we cannot do a
side-by-side comparison against published Cmax/AUC values; the table below is
the simulated NCA summary at the bundle's regimen.

```{r pknca}
sim_nca <- sim |>
  dplyr::filter(time <= sim_hours, !is.na(Cc_col)) |>
  dplyr::transmute(id = id, time = time, Cc_col = Cc_col,
                   treatment = "9MIU_load_then_3MIU_q8h")

dose_df <- cohort |>
  dplyr::filter(evid == 1) |>
  dplyr::transmute(id = id, time = time, amt = amt,
                   treatment = "9MIU_load_then_3MIU_q8h")

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

# At-steady-state interval over the final 24 h of the regimen
intervals <- data.frame(
  start       = sim_hours - 24,
  end         = sim_hours,
  cmax        = TRUE,
  cmin        = TRUE,
  tmax        = TRUE,
  auclast     = TRUE
)

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

The PKNCA output (per-subject Cmax / Cmin / AUClast over the final 24-h
interval) provides the reusable steady-state exposure summary.

# Assumptions and deviations

* **Source publication not on disk.** The Leuppi-Taegtmeyer 2019 AAC
  publication was not present in the operator's literature directory at
  extraction time. All parameter values and equations were taken from the
  DDMORE bundle's `Output_real_*.lst` (final estimates) and
  `Executable_*.mod`. The bundle's `Model_Accommodations.txt` states "There
  are no model differences with the publication referenced", so the published
  final model is faithfully reproduced.
* **Hematocrit fallback.** The `.mod` treats `HT == 0` as a missing-value
  sentinel and substitutes the population-typical hematocrit fraction 0.25.
  The model file reproduces this behaviour verbatim
  (`if (HCT_use <= 0) HCT_use <- 0.25`). The bundle's simulated CSV reports
  `HT = .` (NONMEM-NA equivalent to 0 under `IGNORE=#`) for every subject, so
  the fallback is the actually-used branch in the bundle.
* **Bioavailability F1 is a unit-conversion factor, not absorption.** CMS and
  colistin are administered intravenously; `f(central) <- 1155.5 / 1749.8` is
  the molar-mass ratio of colistin to colistimethate sodium (~0.66) used to
  convert the administered dose (mg of CMS Na) into modelled compartment mass
  (mg colistin base equivalents) so that all concentrations report on the
  conventional colistin-base scale. This is the same construction the source
  uses and is documented inline in the model file.
* **Combined residual error shared across CMS and colistin observations.** The
  source `.mod` `$ERROR` block applies a single combined error model
  (`W = sqrt(THETA(7)^2 * IPRED^2 + THETA(8)^2)`, `Y = IPRED + W * EPS(1)`) to
  both observation compartments without case-splitting on `CMT`. The model
  file declares per-output residual-error parameters
  (`propSd` / `addSd` for `Cc`, `propSd_col` / `addSd_col` for `Cc_col`) per
  the nlmixr2lib multi-output convention, all initialised to the same source
  values (0.222 / 0.459). This is a naming-only deviation; behaviour is
  identical to the source.
* **Cumulative-AUC NONMEM compartments dropped.** `.mod` $DES lines 95-96 add
  bookkeeping integrators `DADT(9) = A(1)/V1` and `DADT(10) = A(2)/V2`; these
  are not needed in nlmixr2lib because the validation vignette uses PKNCA for
  AUC. Mass and concentration ODEs are unchanged.
* **Non-canonical CRRT-specific covariates and compartments.** This is a CRRT
  pharmacokinetic model. The CRRT operational covariates `QBL` (blood flow,
  mL/min) and `QEFF` (effluent flow, mL/h) are not in
  `inst/references/covariate-columns.md` (the canonical register); they are
  declared in this model's `covariateData` and produce
  `checkModelConventions()` warnings of category `covariates`. Similarly the
  CRRT mechanism states `filter`, `cart`, `effl` (and their `_col` colistin
  counterparts) are not canonical compartment names and produce
  `compartments` warnings; they correspond to the CRRT filter / cartridge /
  effluent reservoir physical compartments specific to this model class. The
  metabolite suffix `_col` for colistin is added to
  `R/conventions.R::registeredMetabolites` in the same PR so that
  `central_col`, `Cc_col`, `propSd_col`, and `addSd_col` are canonical.
* **Hematocrit covariate units.** The canonical `HCT` register entry
  documents units in % (0-100). The `.mod` here uses HCT as a fraction (0-1),
  so the model file works with HCT in fractional form. When a user-supplied
  data column carries hematocrit in %, divide by 100 before passing into the
  event table OR rescale inside `model()` before commit-time. This is
  recorded in `covariateData[[HCT]]$notes`.
* **Detailed baseline demographics not reproduced.** Without the publication,
  population fields beyond what the DDMORE bundle's simulated CSV reveals
  (n = 10, adults on CRRT) cannot be filled in faithfully; if the paper PDF
  becomes available the population metadata should be enriched.
