---
title: "Decitabine (Han 2015)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Decitabine (Han 2015)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

mod_meta <- rxode2::rxode(readModelDb("Han_2015_decitabine"))
```

# Model and source

* Citation: `r mod_meta$meta$reference`
* Description: `r mod_meta$meta$description`
* Article (open access): <https://doi.org/10.1186/s13045-015-0208-3>
* PMID: 26482429
* ClinicalTrials.gov: NCT01277484

Decitabine is a hypomethylating cytidine analogue. Han et al. (2015) report
an adaptive phase I trial of post-allogeneic-hematopoietic-stem-cell-
transplantation (allo-HSCT) decitabine maintenance in 16 Korean adults
with higher-risk myelodysplastic syndrome (MDS) or secondary acute myeloid
leukemia (sAML) evolving from MDS. The PK-PD model was used to
individualize the Cycle 2-4 doses for each patient and to estimate the
initial doses for cohorts 2-5. The popPK-PD analysis excluded one female
patient with immune thrombocytopenia (managed with steroids), leaving 15
patients in the mixed-effect analysis (95 PK observations + 311 ANC + 311
PC observations).

The structural model has three layers (Han 2015 Results, "Overall
mixed-effect PK-PD analysis"):

- **PK** -- two-compartment IV-infusion model parameterised on a
  body-surface-area-normalized (per m^2) scale: CL (L/h per m^2),
  Vc (L per m^2), Vp (L per m^2), Q (L/h per m^2). The trial
  administered decitabine as a 60-minute IV infusion at 4-12 mg/m^2/day
  for five consecutive days, repeated every ~4 weeks (mean actual
  interval 34.5 days, SD 8.7).
- **Absolute neutrophil count (ANC)** -- Friberg / Wallin transit
  compartment model with feedback (1 proliferating + 3 transit + 1
  circulating; uniform ktr_N applied across proliferation, transit, and
  circulating-pool clearance). Linear drug effect on the proliferation
  rate (1 - SLOPE_N * Cc) and feedback (BASE_N / circ_anc)^GAMMA_N.
  Baseline is time-invariant: the paper reports "baseline cell count
  increase was not meaningful" for neutrophils.
- **Platelet count (PC)** -- Friberg / Wallin transit compartment model
  with feedback, identical structure to the ANC chain but with a
  time-dependent baseline that increases asymptotically over treatment
  cycles: BASE_P_t = BASE_P + IMP * (1 - exp(-IMK * TIME)). The
  time-dependent BASE_P_t substitutes for BASE_P inside the feedback
  term only; the transit kinetics use the structural ktr_P unchanged.

# Population

```{r population}
mod_meta$meta$population
```

# Source trace

The per-parameter origin is recorded as an in-file comment next to each
`ini()` entry. Every numeric value below was taken from Table 3 of Han
et al. 2015 ("Final parameter estimates and bootstrap outcomes"). The
BSV columns (reported as %CV) were converted to log-scale variances via
omega^2 = log(1 + CV^2); residual error variances were converted to SDs
via sqrt(.) with a /mm^3 -> 10^9/L unit conversion (divide by 1000).

| Equation / parameter | Value | Source location (Han 2015) |
|----------------------|-------|-----------------------------|
| 2-cmt IV infusion ODE  | n/a   | Methods + Results "A two-compartment model was found to best describe the PK data" |
| Friberg / Wallin chain | n/a   | Results, formula `dA(N)/dt = ktr*A(N-1) - ktr*A(N)` + feedback; cites Wallin 2009 [33] |
| `lcl = log(87.8)` (L/h per m^2)   | 87.8     | Table 3, PK parameters row CL |
| `lvc = log(18.5)` (L per m^2)     | 18.5     | Table 3, PK parameters row Vc |
| `lvp = log(22.9)` (L per m^2)     | 22.9     | Table 3, PK parameters row Vp |
| `lq  = log(13.1)` (L/h per m^2)   | 13.1     | Table 3, PK parameters row Q  |
| `lktr_anc   = log(0.0132)` (1/h)  | 0.0132   | Table 3, neutrophil row ktr,N |
| `lslope_anc = log(0.263)` (L/mg)  | 0.263    | Table 3, neutrophil row SLOPE_N |
| `lbase_anc  = log(3.24)` (10^9/L) | 3240/mm^3| Table 3, neutrophil row BASE_N |
| `lgamma_anc = log(0.193)`         | 0.193    | Table 3, neutrophil row GAMMA_N |
| `lktr_plt   = log(0.0244)` (1/h)  | 0.0244   | Table 3, platelet row ktr,P |
| `lslope_plt = log(0.0656)` (L/mg) | 0.0656   | Table 3, platelet row SLOPE_P |
| `lbase_plt  = log(49.2)` (10^9/L) |49200/mm^3| Table 3, platelet row BASE_P |
| `lgamma_plt = log(0.304)`         | 0.304    | Table 3, platelet row GAMMA_P |
| `limp       = log(55)` (10^9/L)   |55000/mm^3| Table 3, platelet row IMP |
| `limk       = log(5.3e-4)` (1/h)  | 0.000530 | Table 3, platelet row IMK |
| `etalcl`        ~ log(1+0.214^2)  | 0.0448   | Table 3, BSV CL = 21.4% CV |
| `etalslope_anc` ~ log(1+0.575^2)  | 0.286    | Table 3, BSV SLOPE_N = 57.5% CV |
| `etalbase_anc`  ~ log(1+0.435^2)  | 0.173    | Table 3, BSV BASE_N = 43.5% CV |
| `etalgamma_anc` ~ log(1+0.394^2)  | 0.144    | Table 3, BSV GAMMA_N = 39.4% CV |
| `etalslope_plt` ~ log(1+0.257^2)  | 0.0640   | Table 3, BSV SLOPE_P = 25.7% CV |
| `etalbase_plt`  ~ log(1+1.050^2)  | 0.743    | Table 3, BSV BASE_P = 105% CV |
| `etalimp`       ~ log(1+0.783^2)  | 0.478    | Table 3, BSV IMP = 78.3% CV |
| `propSd    = sqrt(0.441)`         | 0.664    | Table 3, sigma^2_PK = 0.441 |
| `addSd_ANC = sqrt(754)/1000`      | 0.0275   | Table 3, sigma^2_PD,N = 754 /mm^3^2; /1000 -> 10^9/L |
| `addSd_PLT = sqrt(25000)/1000`    | 0.158    | Table 3, sigma^2_PD,P = 25000 /mm^3^2; /1000 -> 10^9/L |

Table 3 reports "NE" (not estimated) BSV for Vc, Vp, Q, ktr_N, ktr_P,
GAMMA_P, and IMK; these are carried as fixed structural parameters in
the model (no corresponding eta). The paper's covariate analysis
(Results, "No meaningful covariate was found in either the patient
demographic or clinical variables") found no retained covariates for
sex, age, baseline body weight, body-surface-area, or any clinical-lab
variable, so the model has `covariateData = list()`.

# Virtual cohort

The trial schedule for the 15 popPK-PD subjects was decitabine
4-12 mg/m^2/day as a 60-minute IV infusion for 5 consecutive days,
repeated every 4 weeks (28-day protocol) up to Cycle 12; the median
maintenance dose was 7 mg/m^2/day. The first reported recommendation is
5 mg/m^2/day for five consecutive days as the appropriate starting dose.

The simulation below builds a 4-cycle (16-week) virtual cohort at the
recommended initial dose of 5 mg/m^2/day for 5 consecutive days per
cycle, with cycles repeated every 28 days (the protocol interval). PK
parameters are on a per-m^2 scale and the simulated doses are entered
in mg/m^2 directly, so the simulated central concentration is in mg/L
(= ug/mL = 1000 ng/mL).

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

n_subjects <- 200L
mg_per_m2  <- 5            # initial dose
dose_h     <- 1            # 60-minute infusion
inf_rate   <- mg_per_m2 / dose_h  # mg/m^2/h
days_per_cycle <- 5
cycle_h       <- 28 * 24   # 4-week protocol interval (hours)
n_cycles      <- 4
total_h       <- n_cycles * cycle_h

# Dense PK sampling on Cycle 1 Day 1 + sparse PD sampling weekly thereafter,
# matching the paper's protocol (PK on Cycle 1 Day 1; PD weekly through
# Cycle 4).
pk_obs_times <- c(0, 1/3, 2/3, 1, 1.5, 2, 3)  # h after first dose
pd_obs_times <- seq(0, total_h, by = 24 * 7)  # weekly through 16 weeks

# Five daily 1-h infusions per cycle, starting at the cycle start time
dose_starts <- as.vector(outer(
  seq(0, by = 24, length.out = days_per_cycle),         # 5 daily doses
  seq(0, by = cycle_h, length.out = n_cycles),          # cycle starts
  FUN = "+"
))

make_events <- function(id) {
  ev_dose <- data.frame(
    id   = id,
    time = dose_starts,
    evid = 1L,
    amt  = mg_per_m2,
    cmt  = "central",
    rate = inf_rate,
    dur  = NA_real_
  )
  ev_obs_Cc  <- data.frame(id = id, time = pk_obs_times,
                           evid = 0L, amt = 0, cmt = "Cc",
                           rate = NA_real_, dur = NA_real_)
  ev_obs_ANC <- data.frame(id = id, time = pd_obs_times,
                           evid = 0L, amt = 0, cmt = "ANC",
                           rate = NA_real_, dur = NA_real_)
  ev_obs_PLT <- data.frame(id = id, time = pd_obs_times,
                           evid = 0L, amt = 0, cmt = "PLT",
                           rate = NA_real_, dur = NA_real_)
  dplyr::bind_rows(ev_dose, ev_obs_Cc, ev_obs_ANC, ev_obs_PLT)
}

events <- dplyr::bind_rows(lapply(seq_len(n_subjects), make_events))
nrow(events); length(unique(events$id))
```

# Simulation

```{r simulate}
mod <- mod_meta

sim <- rxode2::rxSolve(
  mod,
  events = events
) |>
  as.data.frame()

mod_typical <- mod |> rxode2::zeroRe()
sim_typical <- rxode2::rxSolve(
  mod_typical,
  events = events |> dplyr::filter(id == 1L)
) |>
  as.data.frame()
```

# Replicate published figures

## Cmax cross-check (Discussion, page 4)

The paper's Discussion states: "the predicted C_max after 1-h infusion
of 5 mg/m^2 was 66.0 ng/mL." Reproduce by reading off the typical-value
trajectory at the end of the first infusion (t = 1 h).

**Paper-internal discrepancy.** The published Table 3 parameters imply
an asymptotic infusion concentration k0/CL = 5/87.8 = 57 ng/mL as a
ceiling (i.e. the highest Cc reachable for an infinite-duration
infusion at this dose rate); the actual 1-h end-of-infusion typical
value computed analytically from the 2-cmt parameters in Table 3 is
~51 ng/mL. The paper's quoted 66 ng/mL therefore cannot come from the
Table 3 typical-value PK; it may reflect an alternative dose or rate
than the text states, an inclusion of IIV (a higher percentile of the
predicted-Cmax distribution), or a transcription error in the paper.
The model file in this package faithfully reproduces Table 3, so the
simulated Cmax below matches the parameters as published (~51 ng/mL).
This is flagged as a > 20% discrepancy with the paper's narrative; see
the verification-checklist policy of "flag, do not tune."

```{r cmax-cross-check}
cmax_check <- sim_typical |>
  dplyr::filter(time >= 0, time <= 3) |>
  dplyr::select(time, Cc) |>
  dplyr::mutate(Cc_ngmL = Cc * 1000)

ggplot(cmax_check, aes(time, Cc_ngmL)) +
  geom_line(linewidth = 0.9) +
  geom_vline(xintercept = 1, linetype = 2) +
  labs(x = "Time after first dose (h)",
       y = "Decitabine plasma concentration (ng/mL)",
       title = "Cycle 1 Day 1 typical-value profile after 1-h infusion of 5 mg/m^2",
       caption = "Vertical dashed line marks end of infusion.") +
  theme_minimal()

cmax_pred <- cmax_check |>
  dplyr::filter(time <= 1) |>
  dplyr::pull(Cc_ngmL) |>
  max()

knitr::kable(data.frame(
  metric   = c("Predicted Cmax (ng/mL) at end of 1-h infusion of 5 mg/m^2"),
  this_sim = round(cmax_pred, 1),
  paper    = "~66.0"
))
```

## Figure 2 -- ANC trajectory at 5 mg/m^2/day maintenance over 4 cycles

Han 2015 Figure 2 displays simulated time courses of ANC under the
maintenance dosage of 5 mg/m^2/day for four treatment cycles, including
the median and prediction intervals. The two qualitative features the
text calls out are (i) the lowest ANC value appears in the second cycle
(6-7 weeks after treatment initiation), and (ii) the time to recovery
to ANC > 1000/mm^3 is approximately 5 weeks for the ANC.

```{r figure-2}
sim_anc <- sim |>
  dplyr::filter(!is.na(ANC)) |>
  dplyr::mutate(week = time / 168)

anc_quants <- sim_anc |>
  dplyr::group_by(week) |>
  dplyr::summarise(
    p10 = stats::quantile(ANC, 0.10, na.rm = TRUE),
    p25 = stats::quantile(ANC, 0.25, na.rm = TRUE),
    p50 = stats::quantile(ANC, 0.50, na.rm = TRUE),
    p75 = stats::quantile(ANC, 0.75, na.rm = TRUE),
    p90 = stats::quantile(ANC, 0.90, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(anc_quants, aes(week, p50)) +
  geom_ribbon(aes(ymin = p10, ymax = p90), fill = "steelblue", alpha = 0.15) +
  geom_ribbon(aes(ymin = p25, ymax = p75), fill = "steelblue", alpha = 0.25) +
  geom_line(linewidth = 0.9, colour = "steelblue") +
  geom_hline(yintercept = 1.0, linetype = 2, colour = "darkred") +   # 1000/mm^3 = 1.0 * 10^9/L
  geom_hline(yintercept = 0.5, linetype = 2, colour = "indianred") + # 500/mm^3 = 0.5 * 10^9/L grade 4 toxicity
  labs(x = "Time since treatment initiation (weeks)",
       y = "Absolute neutrophil count (10^9 cells/L)",
       title = "Replicates Han 2015 Figure 2 (ANC time course at 5 mg/m^2/day)",
       caption = paste0(
         "Median, 50% (dark) and 80% (light) prediction intervals from ",
         n_subjects, " subjects. ",
         "Dashed lines: ANC = 1.0 (recovery threshold) and ANC = 0.5 ",
         "(grade 4 toxicity) x 10^9/L."
       )) +
  theme_minimal()

# Time to nadir
nadir_week <- anc_quants$week[which.min(anc_quants$p50)]
recovery_week <- anc_quants |>
  dplyr::filter(week > nadir_week, p50 >= 1.0) |>
  dplyr::slice_head(n = 1) |>
  dplyr::pull(week)

knitr::kable(data.frame(
  metric        = c("Week of median-ANC nadir (Cycle 2 region per paper text)",
                    "Week of median ANC recovery to >= 1.0 x 10^9/L"),
  this_sim      = c(round(nadir_week, 1),
                    if (length(recovery_week)) round(recovery_week, 1) else NA),
  paper_text    = c("3.5 weeks (time-to-nadir), nadir achieved in 2nd cycle (6-7 weeks)",
                    "~5 weeks from last dose")
))
```

## Platelet trajectory with cycle-asymptotic baseline rise

The paper's PD model adds an asymptotic baseline-rise structure for
platelets (BASE_P_t = BASE_P + IMP * (1 - exp(-IMK * TIME))) to capture
the gradual platelet-count increase over cycles -- a contribution of
decitabine to cell proliferation reported in previous studies.

```{r figure-platelet}
sim_plt <- sim |>
  dplyr::filter(!is.na(PLT)) |>
  dplyr::mutate(week = time / 168)

plt_quants <- sim_plt |>
  dplyr::group_by(week) |>
  dplyr::summarise(
    p10 = stats::quantile(PLT, 0.10, na.rm = TRUE),
    p25 = stats::quantile(PLT, 0.25, na.rm = TRUE),
    p50 = stats::quantile(PLT, 0.50, na.rm = TRUE),
    p75 = stats::quantile(PLT, 0.75, na.rm = TRUE),
    p90 = stats::quantile(PLT, 0.90, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(plt_quants, aes(week, p50)) +
  geom_ribbon(aes(ymin = p10, ymax = p90), fill = "darkgreen", alpha = 0.15) +
  geom_ribbon(aes(ymin = p25, ymax = p75), fill = "darkgreen", alpha = 0.25) +
  geom_line(linewidth = 0.9, colour = "darkgreen") +
  geom_hline(yintercept = 30, linetype = 2, colour = "darkred") +  # 30,000/mm^3 holding threshold
  geom_hline(yintercept = 25, linetype = 3, colour = "indianred") +# 25,000/mm^3 grade 4 toxicity
  labs(x = "Time since treatment initiation (weeks)",
       y = "Platelet count (10^9 cells/L)",
       title = "Platelet time course at 5 mg/m^2/day maintenance",
       caption = paste(
         "Median + 50% (dark) and 80% (light) prediction intervals.",
         "Dashed line: PLT = 30 x 10^9/L (dose-hold);",
         "dotted line: PLT = 25 x 10^9/L (grade 4 toxicity)."
       )) +
  theme_minimal()
```

# PKNCA validation

The paper's Discussion provides two cross-check NCA values: predicted
Cmax = 66.0 ng/mL after a 1-h infusion of 5 mg/m^2, and a typical
half-life of 0.31 h. PKNCA below evaluates Cycle-1-Day-1 (single-dose)
NCA on the dense PK sampling grid (0, 20, 40, 60, 90, 120, 180 min)
and compares Cmax and half-life with the paper.

```{r pknca}
sim_nca <- sim |>
  dplyr::filter(time <= 3, !is.na(Cc)) |>
  dplyr::distinct(id, time, .keep_all = TRUE) |>
  dplyr::select(id, time, Cc) |>
  dplyr::mutate(Cc_ngmL = Cc * 1000)

dose_df <- events |>
  dplyr::filter(evid == 1L, time == 0) |>
  dplyr::distinct(id, time, .keep_all = TRUE) |>
  dplyr::transmute(id, time, amt, treatment = "5 mg/m^2 1-h IV")

conc_df <- sim_nca |> dplyr::mutate(treatment = "5 mg/m^2 1-h IV")

conc_obj <- PKNCA::PKNCAconc(
  conc_df,
  Cc_ngmL ~ time | treatment / id,
  concu = "ng/mL",
  timeu = "h"
)
dose_obj <- PKNCA::PKNCAdose(
  dose_df,
  amt ~ time | treatment + id,
  doseu = "mg/m^2"
)

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

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

knitr::kable(summary(nca_res),
             caption = "Cycle-1 Day-1 single-dose NCA (200 subjects, 5 mg/m^2 1-h IV).")
```

```{r pknca-vs-paper}
pp_res <- as.data.frame(nca_res$result)
cmax_median <- median(pp_res$PPORRES[pp_res$PPTESTCD == "cmax"], na.rm = TRUE)
hl_median   <- median(pp_res$PPORRES[pp_res$PPTESTCD == "half.life"], na.rm = TRUE)

knitr::kable(data.frame(
  metric        = c("Cmax (ng/mL)", "Apparent half-life (h)"),
  this_sim_med  = c(round(cmax_median, 1), round(hl_median, 2)),
  paper_target  = c("~66 (Discussion, page 4)", "~0.31 (Discussion, page 4)")
))
```

# Assumptions and deviations

- **Per-m^2 PK scale and dosing unit.** Han 2015 reports PK parameters
  on a body-surface-area-normalized scale (L/h per m^2 and L per m^2);
  this is the convention for cytotoxic dosing in oncology. Doses are
  consequently entered in `mg/m^2` (not `mg`), and the central
  compartment carries `mg/m^2` -- `central / vc` yields `mg/L` directly.
  Users who wish to dose in `mg` for a patient with known BSA must
  pre-multiply the absolute mg dose by `1 / BSA` (m^2) before passing
  it to the model, or alternatively rescale the per-m^2 PK parameters
  by their patient's BSA. The trial protocol always recorded doses in
  `mg/m^2/day`, so the per-m^2 model is the operational convention.
- **No covariates in the final model.** Han 2015 reports "no meaningful
  covariate was found in either the patient demographic or clinical
  variables", so `covariateData = list()`. Demographic facts (age range
  19-64 years, 40% female, all Korean) are recorded in `population`
  metadata for downstream filtering but do not affect simulation.
- **Vc, Vp, Q, ktr_N, ktr_P, GAMMA_P, and IMK have no IIV.** Table 3's
  BSV column reports "NE" (not estimated) for these parameters; the
  Methods says "BSV ... was the only random effect which could be
  estimated, except for the proportional residual error" for PK CL.
  The model carries these as typical-value-only parameters (no
  corresponding eta in `ini()`), consistent with the paper's reporting.
- **Friberg / Wallin transit kinetics use a single ktr.** Both ANC and
  PLT chains use one structural ktr (the rate constant of inter-
  compartmental movement, reported in 1/h directly in Table 3) applied
  uniformly to the proliferation, transit, and circulating
  compartments. This matches the Wallin 2009 form the paper cites
  (reference 33). The implied mean transit time MTT = 4 / ktr is
  303 h (12.6 days) for ANC and 164 h (6.8 days) for PLT.
- **Time-varying platelet baseline enters via feedback only.** The
  `BASE_P_t = BASE_P + IMP * (1 - exp(-IMK * TIME))` structure
  substitutes for BASE_P inside the Friberg feedback term
  `(BASE_P_t / circ_plt)^GAMMA_P` only; the transit kinetics use the
  structural `ktr_P` unchanged, and the initial conditions
  `precursor*_plt(0) = circ_plt(0) = BASE_P` are correct because
  `BASE_P_t(t = 0) = BASE_P`. This faithfully reproduces the paper's
  second formula in "Overall mixed-effect PK-PD analysis".
- **Time origin for IMK.** The paper writes the platelet asymptote in
  terms of `TIME` (time from initiation of decitabine treatment). The
  model uses rxode2's reserved variable `t`, which is simulation time;
  for any rxode2 events table that places the first dose at `t = 0`
  the two are identical. Users who define event tables with t < 0 for
  pre-dose observations will see `BASE_P_t < BASE_P` because
  `1 - exp(-IMK*t) < 0` for `t < 0`; the canonical usage is to start
  at `t = 0`.
- **Initial-condition simplification.** Friberg chains are placed at
  steady-state at `t = 0` (all five compartments equal to BASE,
  consistent with the Friberg derivation where the unperturbed system
  has `precursor1 = ... = circ` at the baseline value). The paper does
  not provide a separate ramp-in period; the trial enrolled patients
  whose blood counts had recovered to PC > 30,000/mm^3 and ANC >
  1000/mm^3 prior to decitabine initiation.
- **Half-life cross-check.** The paper's Discussion reports an "average
  terminal half-life ... of 0.31 h in this study". With CL/Vc = 87.8 /
  18.5 = 4.745 /h the rapid distribution half-life is `ln(2)/4.745 =
  0.146 h`; the slower beta phase is around 1.4 h (per analytic 2-cmt
  eigenvalues from CL, Vc, Vp, Q). The simulated single-dose half-life
  PKNCA estimates (above) reflect the average behaviour across the
  observed sampling window (3 h post-dose) where the alpha and beta
  phases mix; the apparent half-life depends sensitively on which
  terminal-phase points the NCA algorithm uses. The published "0.31 h"
  is an empirical estimate from non-compartmental analysis on the
  actual clinical samples (Discussion, p. 4) -- the PKNCA value here is
  on simulated typical-value data and is reported only as a sanity
  check, not as a faithful reproduction of the empirical estimate.
- **`circ_anc`, `precursor*_anc`, `circ_plt`, `precursor*_plt`
  compartment names deviate from the canonical compartment register**
  (`depot`, `central`, `peripheral<n>`, `effect`, ...). They are
  retained because they directly mirror the standard Friberg / Wallin
  notation (proliferation -> 3 transits -> circulating); the
  Friberg-precedent files (`Friberg_2002_paclitaxel` in DDMORE,
  `Petrov_2024_romiplostim`, `Netterberg_2017_docetaxel`,
  `Guo_2022_PF_06939999`) use the same naming.
  `checkModelConventions()` flags these as warnings of
  `compartments` severity; the deviation is intentional.
- **Residual error stored in `ini()` on the canonical 10^9/L PD unit.**
  Han 2015 reports the PD additive residual variances as
  `sigma^2_PD,N = 754 /mm^3^2` and `sigma^2_PD,P = 25,000 /mm^3^2`. The
  SDs in `mm^-3` are `sqrt(754) = 27.46` and `sqrt(25000) = 158.1`. To
  match `addSd_ANC` / `addSd_PLT` to the canonical PD observation unit
  of `10^9 cells/L` (i.e. `1000 / mm^3`), the SDs in `ini()` are
  divided by 1000 to give `0.0275` and `0.158`.
- **Successful-bootstrap convergence not propagated.** Han 2015 reports
  bootstrap convergence rates of 78.8% (PK) and 78.0% (PD) with median
  parameter estimates that closely match the point estimates (Table 3,
  Bootstrap median column). The model carries the point estimates
  rather than the bootstrap medians; the differences are < 5% for
  every reported parameter.
