---
title: "MTMSA-Trp (Niloy 2026)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{MTMSA-Trp (Niloy 2026)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: `r rxode2::rxode(readModelDb("Niloy_2026_MTMSATrp_mouse"))$reference`
* Description: `r rxode2::rxode(readModelDb("Niloy_2026_MTMSATrp_mouse"))$description`
* Article: <https://doi.org/10.21203/rs.3.rs-9035594/v1>

MTMSA-Trp is a novel mithramycin analogue under investigation for Ewing
sarcoma. Niloy 2026 develops a one-compartment population PK model from
sparse and destructive plasma sampling in female athymic nu/nu mice
following single IV bolus administration over a 30-fold dose range
(0.3 to 10 mg/kg). Apparent clearance decreases with increasing dose;
the authors describe this as an empirical (not mechanistic) feature.

# Population

The MTMSA-Trp popPK dataset comprised 121 plasma concentrations from 70
female athymic nu/nu mice (Envigo RMS) following single IV bolus
administration of 0.3, 1, 3, 5, or 10 mg/kg. Sampling combined serial
and terminal / destructive designs, so some animals contributed a single
concentration. Plasma was quantified by a validated LC-MS/MS bioanalytical
method. Mice were healthy and xenograft-naive (the PK dataset is
distinct from the efficacy / xenograft experiments). All work was
performed under approved IACUC protocols at St. Jude Children's Research
Hospital and the University of Kentucky (Niloy 2026 Methods, "Mice" and
"Pharmacokinetic data").

The same information is available programmatically via the model's
`population` metadata
(`readModelDb("Niloy_2026_MTMSATrp_mouse")$population` after the model
is loaded).

# Source trace

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

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lcl <- log(CL)`     | CL = 39.18 mL/h/kg at 3 mg/kg | Niloy 2026 Table 1 (CL, RSE 4.97%) |
| `lvc <- log(V)`      | V  = 53.06 mL/kg              | Niloy 2026 Table 1 (V,  RSE 7.55%) |
| `e_cohdose_cl` (beta) | -0.30                        | Niloy 2026 Table 1 (beta, RSE 10.9%) |
| `etalcl` (omega^2)   | log(1 + 0.058^2) = 0.003359  | Niloy 2026 Table 1 (IIV on CL = 5.8% CV) |
| `etalvc` (omega^2)   | log(1 + 0.25^2)  = 0.060625  | Niloy 2026 Table 1 (IIV on V  = 25% CV) |
| `propSd`             | 0.51 (51% proportional CV)   | Niloy 2026 Table 1 (b, RSE 8.13%) |
| Reference dose       | 3 mg/kg                      | Niloy 2026 Methods, "Model development", Eq. 2 |
| Lognormal IIV form   | `P_i = P_pop * exp(eta_i)`   | Niloy 2026 Methods (exponential IIV model) |
| Power covariate form | `CL_i = CL_pop * (Dose_i / 3)^beta` | Niloy 2026 Methods, Eq. 2 |
| Residual error form  | proportional                 | Niloy 2026 Methods / Results (selected over additive and combined) |
| Structural model     | 1-compartment, first-order elim. | Niloy 2026 Results, "Model development" |

# Virtual cohort

The original observed concentrations are not publicly available. The
figures below use a virtual cohort with the same five dose levels and
study design as Niloy 2026.

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

# Sampling grid: dense early (0-2 h) to characterize the C0 drop on an
# IV bolus with t1/2 ~ 0.5-1.5 h, then sparser through 24 h to
# capture the terminal phase fully (~16-50 half-lives across the dose
# range -- enough for AUC_0-inf).
obs_times <- c(0, 0.083, 0.25, 0.5, 1, 2, 4, 6, 8, 12, 16, 24)

doses <- c(0.3, 1, 3, 5, 10)  # mg/kg, the five Niloy 2026 cohorts
n_per_dose <- 200L

make_cohort <- function(dose_mgkg, n, id_offset) {
  ids <- id_offset + seq_len(n)
  dose_rows <- tibble::tibble(
    id      = ids,
    time    = 0,
    amt     = dose_mgkg,    # per-kg dose carried directly on the event record
    evid    = 1L,
    cmt     = "central",
    COHDOSE = dose_mgkg
  )
  obs_rows <- tidyr::expand_grid(
    id      = ids,
    time    = obs_times
  ) |>
    dplyr::mutate(
      amt     = NA_real_,
      evid    = 0L,
      cmt     = "central",  # ODE-state name (not the algebraic observable "Cc")
      COHDOSE = dose_mgkg
    )
  dplyr::bind_rows(dose_rows, obs_rows) |>
    dplyr::arrange(id, time, dplyr::desc(evid)) |>
    dplyr::mutate(treatment = sprintf("%.1f mg/kg", dose_mgkg))
}

events <- purrr::map_dfr(
  seq_along(doses),
  function(i) make_cohort(doses[i], n_per_dose,
                          id_offset = (i - 1L) * n_per_dose)
)

# Disjoint-ID assertion -- guards against the silent multi-cohort
# Frankenstein-subject bug in rxSolve (see vignette-template Notes).
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
```

# Simulation

```{r simulate}
mod <- readModelDb("Niloy_2026_MTMSATrp_mouse")

sim <- rxode2::rxSolve(mod,
                       events = events,
                       keep   = c("treatment", "COHDOSE")) |>
  as.data.frame()
```

For deterministic replication (reproducing the typical-value profile
without between-subject variability), zero out the random effects:

```{r simulate-typical}
mod_typical <- mod |> rxode2::zeroRe()
sim_typical <- rxode2::rxSolve(mod_typical,
                               events = events,
                               keep   = c("treatment", "COHDOSE")) |>
  as.data.frame()
```

# Replicate published figures

## Figure 1 -- observed concentration-time profiles by dose

```{r figure-1}
# Replicates Figure 1 of Niloy 2026: median PK profile by dose group.
sim_typical |>
  dplyr::filter(id == 1 + (match(round(COHDOSE, 3), round(doses, 3)) - 1L) * n_per_dose,
                time > 0) |>
  ggplot(aes(time, Cc, colour = treatment)) +
  geom_line(linewidth = 0.7) +
  scale_y_log10() +
  labs(x = "Time (h)", y = "MTMSA-Trp plasma (ng/mL)",
       colour = "Dose",
       title = "Figure 1 -- typical-value PK profile by dose",
       caption = "Replicates Figure 1 of Niloy 2026 (mean profile across the dose range).")
```

## Figure 3 -- VPC of the final model

```{r figure-3}
# Replicates Figure 3 of Niloy 2026: pooled VPC with 5th / 50th / 95th
# percentiles over time.
sim |>
  dplyr::filter(time > 0) |>
  dplyr::group_by(time, treatment) |>
  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.7) +
  facet_wrap(~ treatment, scales = "free_y") +
  scale_y_log10() +
  labs(x = "Time (h)", y = "MTMSA-Trp plasma (ng/mL)",
       title = "Figure 3 -- VPC of final popPK model",
       caption = paste("Replicates Figure 3 of Niloy 2026; ribbon = 5th-95th,",
                       "line = median across", n_per_dose, "virtual mice per dose."))
```

## Figure 5 -- dose-dependent nonlinearity

```{r figure-5}
# Replicates Figure 5 of Niloy 2026: simulated AUC_0-inf and
# dose-normalised AUC_0-inf versus dose. Use the typical-value
# simulation (rxode2::zeroRe) for the population-mean exposure across
# the dose range. AUC_0-inf for a 1-compartment IV bolus simplifies to
#   AUC_inf = Dose / CL
# evaluated at the individual's CL. For a typical mouse this is
# analytic, but here we integrate the simulated profile via PKNCA-style
# trapezoidal AUC and extrapolation, then plot the median +/- 5-95
# percentile band from the stochastic simulation as paper does.
auc_per_subject <- sim |>
  dplyr::filter(time > 0) |>
  dplyr::group_by(id, treatment, COHDOSE) |>
  dplyr::arrange(time) |>
  dplyr::summarise(
    auc_last = sum(diff(time) * (head(Cc, -1) + tail(Cc, -1)) / 2),
    clast    = dplyr::last(Cc),
    tlast    = dplyr::last(time),
    # Use the model's typical-value kel for extrapolation: kel = CL / V
    # at the dose-adjusted CL. Per-subject extrapolation is robust here
    # because the model is structural and noise-free apart from IIV.
    kel_typ  = (39.18 * (dplyr::first(COHDOSE) / 3)^(-0.30)) / 53.06,
    auc_inf  = auc_last + clast / kel_typ,
    .groups  = "drop"
  )

auc_summary <- auc_per_subject |>
  dplyr::group_by(treatment, COHDOSE) |>
  dplyr::summarise(
    auc_med = median(auc_inf),
    auc_lo  = quantile(auc_inf, 0.05),
    auc_hi  = quantile(auc_inf, 0.95),
    .groups = "drop"
  ) |>
  dplyr::mutate(
    dn_med = auc_med / COHDOSE,
    dn_lo  = auc_lo  / COHDOSE,
    dn_hi  = auc_hi  / COHDOSE
  )

p_auc <- ggplot(auc_summary, aes(COHDOSE, auc_med)) +
  geom_errorbar(aes(ymin = auc_lo, ymax = auc_hi), width = 0.1) +
  geom_point(size = 2) +
  scale_x_log10() +
  scale_y_log10() +
  labs(x = "Dose (mg/kg)", y = "AUC_0-inf (ng*h/mL)",
       title = "Figure 5a -- exposure vs dose")

p_dn  <- ggplot(auc_summary, aes(COHDOSE, dn_med)) +
  geom_errorbar(aes(ymin = dn_lo, ymax = dn_hi), width = 0.1) +
  geom_point(size = 2) +
  scale_x_log10() +
  labs(x = "Dose (mg/kg)", y = "Dose-normalised AUC_0-inf (ng*h/mL per mg/kg)",
       title = "Figure 5b -- dose-normalised exposure vs dose")

p_auc; p_dn
```

# PKNCA validation

```{r pknca}
# Keep only observation rows; the filter is intentionally minimal
# (only !is.na(Cc)) so the time=0 row is preserved for AUC anchoring.
sim_nca <- sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, treatment)

# Defensive time=0 guarantee: an IV bolus produces a non-zero C(0+);
# PKNCA wants a time=0 row to anchor AUC. If the simulation grid
# already includes time = 0, distinct() collapses duplicates.
sim_nca <- dplyr::bind_rows(
  sim_nca,
  sim_nca |> dplyr::distinct(id, treatment) |>
    dplyr::mutate(time = 0,
                  Cc   = first(sim_nca$Cc[sim_nca$time > 0]) * 0 + NA_real_)
) |>
  dplyr::group_by(id, treatment, time) |>
  dplyr::slice_head(n = 1) |>
  dplyr::ungroup() |>
  dplyr::arrange(id, treatment, time)

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | treatment + id,
                             concu = "ng/mL", timeu = "h")

dose_df <- events |>
  dplyr::filter(evid == 1L) |>
  dplyr::select(id, time, amt, treatment)

dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id,
                             doseu = "mg/kg")

intervals <- data.frame(
  start       = 0,
  end         = Inf,
  cmax        = TRUE,
  tmax        = TRUE,
  aucinf.obs  = TRUE,
  half.life   = TRUE
)

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

## Comparison against published / analytic NCA

Niloy 2026 reports only typical PK parameters (CL, V, beta) and does
not tabulate per-dose NCA values; the dose-stratified results appear
graphically in Figure 5. The "Reference" column below is the analytic
expectation derived from Niloy 2026 Table 1 parameters:

- Cmax = Dose / V    (mg/mL converted to ng/mL = x 1e6)
- AUC_0-inf = Dose / CL_i, with CL_i = CL_pop * (Dose / 3)^beta
- Half-life = ln(2) * V / CL_i

```{r nca_compare}
analytic_reference <- function(dose_mgkg) {
  CL_pop <- 39.18; V <- 53.06; beta <- -0.30
  CL_i   <- CL_pop * (dose_mgkg / 3)^beta
  tibble::tibble(
    treatment   = sprintf("%.1f mg/kg", dose_mgkg),
    cmax        = (dose_mgkg / V) * 1e6,
    tmax        = 0,
    aucinf.obs  = (dose_mgkg / CL_i) * 1e6,
    half.life   = log(2) * V / CL_i
  )
}

published <- purrr::map_dfr(doses, analytic_reference)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_res,
  reference = published,
  by        = "treatment",
  units     = c(cmax       = "ng/mL",
                aucinf.obs = "ng*h/mL",
                tmax       = "h",
                half.life  = "h"),
  tolerance_pct = 20
)

knitr::kable(
  cmp,
  caption = paste("Simulated (median across", n_per_dose,
                  "virtual mice) vs. analytic prediction from Niloy 2026",
                  "Table 1. * differs by >20%."),
  align   = c("l", "l", "r", "r", "r")
)
```

Any starred row should be investigated but not tuned -- they typically
reflect the high IIV / shrinkage on CL acknowledged in Niloy 2026
(Discussion).

# Assumptions and deviations

- **Per-kg parameterisation kept intact.** Niloy 2026 reports CL in
  mL/h/kg and V in mL/kg with dose in mg/kg. The model file leaves the
  parameters in these per-kg units; the event-record `amt` carries the
  per-kg dose directly. There is no explicit body-weight covariate,
  consistent with the source's choice not to model weight as a
  covariate (no individual mouse weights are tabulated).
- **Concentration unit set to ng/mL.** The native units of the model
  (mg/kg dose / mL/kg V) give Cc in mg/mL. The observation expression
  multiplies by 1e6 to express plasma concentration in ng/mL, matching
  the typical LC-MS/MS reporting unit and the implicit scale of
  Niloy 2026 Figure 1. The reference bioanalytical paper cited as [8]
  in Niloy 2026 is not on disk; the units assumption is documented
  here in case it differs.
- **IIV variance encoding.** Niloy 2026 Table 1 reports IIV as %CV
  (5.8% on CL, 25% on V). These were converted to NONMEM-style omega^2
  via the Monolix log-normal convention `omega^2 = log(1 + CV^2)`,
  giving 0.003359 (CL) and 0.060625 (V). The bootstrap median column
  reports slightly different point estimates (0.061 and 0.21 for
  CL / V respectively, on the omega-SD scale); the model uses the
  Estimate column (the final MLE), not the bootstrap median.
- **Residual error coefficient interpretation.** Table 1 reports
  `b = 0.51` under a row labelled `b (%)`. The "(%)" tag appears to be
  carried over from the IIV %CV rows above and does not apply
  literally to the proportional residual coefficient b, which is a
  decimal fraction in Monolix output. The decimal interpretation
  (51% proportional CV) is confirmed by the bootstrap median (0.5) and
  95% CI (0.41-0.58), and is what the model uses.
- **Reference NCA values are analytic, not measured.** The
  Comparison section uses analytic predictions from the Niloy 2026
  parameter estimates (Cmax = Dose/V, AUC_inf = Dose/CL_i,
  t1/2 = ln(2) * V / CL_i) because the paper does not tabulate
  per-dose NCA values; dose-stratified exposure appears only
  graphically (Figure 5). A starred row indicates that the
  simulation median departs from the analytic prediction by more
  than 20%, which is expected at the lower dose levels given the
  reported shrinkage on CL (81%) and is documented in the source
  Discussion rather than something to be tuned.
