---
title: "Oxaliplatin (Valenzuela 2011)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Oxaliplatin (Valenzuela 2011)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: Valenzuela B, Nalda-Molina R, Bretcha-Boix P, Escudero-Ortiz V,
  Duart MJ, Carbonell V, Sureda M, Rebollo JP, Farre J, Brugarolas A,
  Perez-Ruixo JJ. *Pharmacokinetic and Pharmacodynamic Analysis of
  Hyperthermic Intraperitoneal Oxaliplatin-Induced Neutropenia in
  Subjects with Peritoneal Carcinomatosis.* AAPS J. 2011 Mar;13(1):72-82.
* DOI: <https://doi.org/10.1208/s12248-010-9249-2> (PMID 21210260).
* Erratum: AAPS J. 2011 Jun;13(2):318 (referenced in the PubMed
  "Erratum in" field of PMID 21210260; the erratum text was not
  available on disk at extraction time -- see Section Errata).
* Friberg backbone: Friberg LE, Henningsson A, Maas H, Nguyen L,
  Karlsson MO. *Model of chemotherapy-induced myelosuppression with
  parameter consistency across drugs.* J Clin Oncol 2002;20(24):4713-4721,
  <https://doi.org/10.1200/JCO.2002.02.140> (PMID 12488418).

```{r load-model}
mod <- rxode2::rxode(readModelDb("Valenzuela_2011_oxaliplatin"))
cat(mod$meta$description, "\n")
```

# Population

The study enrolled 30 adults with peritoneal carcinomatosis at USP Hospital
San Jaime (Torrevieja, Spain) between 2006 and 2009. Baseline characteristics
are summarized in Table I of the source: mean age 57.9 years (SD 10.5, range
32-75); mean body weight 69.3 kg (SD 12.1, range 42-90); mean BSA 1.7 m^2
(SD 0.2, range 1.4-2.0); 60% female. Primary tumour types were ovarian (n=10),
colorectal (n=9), appendiceal (n=5), gastric (n=3), endometrial (n=2), and
primary papillary (n=1). All patients underwent cytoreductive surgery
followed by hyperthermic intraperitoneal oxaliplatin (HIO) 360 mg/m^2 in
4% icodextrin perfusate (2.5-6.0 L) for 30-60 min at perfusate
temperature 42-43 degC, then intraperitoneal 5-FU 15 mg/kg over a 1-h
infusion on each of the 5 postoperative days. The model uses 140 peritoneal
oxaliplatin concentrations, 338 plasma concentrations, and 678 ANC values.

The same information is available programmatically via
`rxode2::rxode(readModelDb("Valenzuela_2011_oxaliplatin"))$meta$population`.

```{r population-metadata}
str(mod$meta$population, max.level = 1L)
```

# Source trace

The per-parameter origin is recorded as an in-file comment next to each
`ini()` entry in
`inst/modeldb/specificDrugs/Valenzuela_2011_oxaliplatin.R`. The table
below collects the population point estimates and ODE-equation references
in one place for review against the source PDF.

| Equation / parameter                                | Value          | Source location                |
|-----------------------------------------------------|----------------|--------------------------------|
| `Qa` (peritoneum-to-plasma clearance)               | 2.70 L/h       | Table II, Pharmacokinetic block |
| `Va` (peritoneum volume, mapped onto `lvd`)         | 8.33 L         | Table II                       |
| `Cl/F` (apparent systemic clearance)                | 1.61 L/h       | Table II                       |
| `Q2/F` (apparent intercompartmental clearance)      | 77.0 L/h       | Table II                       |
| `Vc/F` (apparent central volume)                    | 19.2 L         | Table II                       |
| `Vp/F` (apparent peripheral volume)                 | 72.8 L         | Table II                       |
| `Circ0` (baseline ANC)                              | 7.05 x 10^9/L  | Table II, System block         |
| `MTT` (mean transit time, FIXED)                    | 118 h          | Table II 'MTT 118 Fixed'; ref 30 (Friberg 2002) |
| `gamma` (feedback exponent)                         | 0.135          | Table II                       |
| `alpha` (linear drug-effect slope)                  | 0.182 L/mg     | Table II, Drug block; Eq 9     |
| `kCirc` (circ elimination rate, FIXED)              | 0.07 1/h       | Methods text + ref 30 (Friberg 2002) |
| `omega_Qa, omega_Va, omega_Cl/F, omega_Vc/F, omega_Vp/F` (CV %) | 34.1, 17.7, 85.6, 57.9, 23.5 | Table II IIV block |
| `omega_Circ0, omega_MTT, omega_alpha` (CV %)        | 42.3, 32.8, 141 | Table II IIV block             |
| `sigma_1` (plasma), `sigma_2` (peritoneum), `sigma` (ANC) (CV %) | 14.7, 16.5, 49.7 | Table II Residual block |
| Eq 1: `dA/dt = -ka * A` (peritoneum, ka = Qa/Va)    | n/a            | Page 75                        |
| Eq 2: `dC/dt = Qa/Va * A - Cl/Vc * C - Q2/Vc * C + Q2/Vp * P` | n/a  | Page 75                        |
| Eq 3: `dP/dt = Q2/Vc * C - Q2/Vp * P`               | n/a            | Page 75                        |
| Eq 4: `dProl/dt = kprol * Prol * (Circ0/Circ)^gamma * (1 - Edrug) - ktr * Prol` | n/a | Page 75 |
| Eq 5-7: `dTransit_i/dt = ktr * Transit_{i-1} - ktr * Transit_i` (i = 1..3) | n/a | Page 75 |
| Eq 8: `dCirc/dt = ktr * Transit3 - kCirc * Circ`    | n/a            | Page 75 (kCirc separate per text) |
| Eq 9: `Edrug = alpha * Cc`                          | n/a            | Page 76                        |
| Eq 10: `MTT = (n+1) / ktr` with n = 3 transits      | n/a            | Page 76                        |

# Virtual cohort

The paper screened ten demographic and biochemistry covariates (age, sex,
body weight, serum creatinine, serum albumin, ALT, AST, total bilirubin,
hemoglobin, hematocrit) and found no correlation with PK / PD parameters,
so the final model carries no covariate effects. The simulation cohort
below uses BSA-derived single-dose amounts drawn from the Table I BSA
distribution; no other subject-level covariates feed into the model
predictions.

```{r cohort}
set.seed(20260603)
n_sub <- 200L

cohort <- tibble::tibble(
  id  = seq_len(n_sub),
  BSA = pmax(1.0, rnorm(n_sub, mean = 1.7, sd = 0.2))   # Table I baseline BSA (m^2)
)
cohort$amt <- 360 * cohort$BSA                          # 360 mg/m^2 single HIO dose

# Build events: bolus dose into peritoneum at t = 0, then dense observation
# grid covering ~28 days (paper's PD follow-up window).
obs_times <- sort(unique(c(
  seq(0, 4, length.out = 60),
  seq(4, 28, by = 2),
  seq(28, 28 * 24, by = 6)
)))

dose_rows <- cohort |>
  dplyr::transmute(
    id, time = 0, amt = amt,
    evid = 1L, cmt = "depot", BSA = BSA
  )

obs_rows <- tidyr::expand_grid(
  id   = cohort$id,
  time = obs_times
) |>
  dplyr::mutate(
    amt = NA_real_, evid = 0L, cmt = "Cc",
    BSA = cohort$BSA[match(id, cohort$id)]
  )

events <- dplyr::bind_rows(dose_rows, obs_rows) |>
  dplyr::arrange(id, time)

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

# Simulation

```{r simulate}
sim <- rxode2::rxSolve(mod, events = events, keep = "BSA",
                       returnType = "data.frame")
sim$day <- sim$time / 24
```

A deterministic typical-value simulation (zeroing the between-subject
random effects) is also useful for replicating the paper's figures.

```{r simulate-typical}
mod_typ <- rxode2::zeroRe(mod)
events_typ <- events |> dplyr::filter(id == 1L)
sim_typ <- rxode2::rxSolve(mod_typ, events = events_typ, keep = "BSA",
                           returnType = "data.frame")
sim_typ$day <- sim_typ$time / 24
```

# Steady-state ANC baseline

Before any drug exposure the proliferating + transit chain holds the
circulating ANC at `Circ0` indefinitely. With `kCirc` distinct from `ktr`
the precursor pools sit at `(kCirc / ktr) * Circ0` (see model file
`prol_ss`); the typical-value baseline should land on the reported
`Circ0 = 7.05 x 10^9/L` and stay there until the HIO bolus drives the
chain.

```{r steady-state}
typ0 <- sim_typ[sim_typ$time == 0, c("time","Cperi","Cc","ANC","precursor1","circ")]
print(typ0)
```

# Replicate Figure 3 (visual predictive check)

Figure 3 of Valenzuela 2011 shows the 5th, 50th, and 95th percentiles of
the simulated peritoneal and plasma concentrations and ANC, with the
observed percentiles overlaid. We reproduce the same percentile bands
from the packaged model below. The observed data themselves are not
public; only the simulation envelope can be shown here.

```{r figure-3-pk, fig.width = 8, fig.height = 5}
pk_vpc <- sim |>
  dplyr::filter(time > 0) |>
  dplyr::group_by(time) |>
  dplyr::summarise(
    Cperi_Q05 = quantile(Cperi, 0.05, na.rm = TRUE),
    Cperi_Q50 = quantile(Cperi, 0.50, na.rm = TRUE),
    Cperi_Q95 = quantile(Cperi, 0.95, na.rm = TRUE),
    Cc_Q05    = quantile(Cc,    0.05, na.rm = TRUE),
    Cc_Q50    = quantile(Cc,    0.50, na.rm = TRUE),
    Cc_Q95    = quantile(Cc,    0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  dplyr::mutate(hours = time)

pk_long <- pk_vpc |>
  tidyr::pivot_longer(
    cols = -c(time, hours),
    names_to = c("compartment", ".value"),
    names_pattern = "(Cperi|Cc)_(Q05|Q50|Q95)"
  ) |>
  dplyr::mutate(
    compartment = dplyr::recode(compartment,
                                Cperi = "Peritoneum",
                                Cc    = "Plasma")
  )

ggplot(pk_long, aes(hours, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25, fill = "steelblue") +
  geom_line(colour = "steelblue") +
  facet_wrap(~compartment, scales = "free_y") +
  scale_y_log10() +
  coord_cartesian(xlim = c(0, 28)) +
  labs(x = "Time (h)", y = "Oxaliplatin concentration (mg/L)",
       title = "Replicates Figure 3 upper / mid panels of Valenzuela 2011",
       caption = "5th / 50th / 95th simulated percentiles, n = 200 virtual subjects.") +
  theme_minimal()
```

```{r figure-3-anc, fig.width = 8, fig.height = 4}
anc_vpc <- sim |>
  dplyr::filter(time > 0) |>
  dplyr::group_by(time) |>
  dplyr::summarise(
    ANC_Q05 = quantile(ANC, 0.05, na.rm = TRUE),
    ANC_Q50 = quantile(ANC, 0.50, na.rm = TRUE),
    ANC_Q95 = quantile(ANC, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  dplyr::mutate(day = time / 24)

ggplot(anc_vpc, aes(day, ANC_Q50)) +
  geom_ribbon(aes(ymin = ANC_Q05, ymax = ANC_Q95), alpha = 0.25, fill = "firebrick") +
  geom_line(colour = "firebrick") +
  geom_hline(yintercept = c(0.5, 1.0), linetype = "dashed", alpha = 0.4) +
  labs(x = "Time (days)", y = "ANC (10^9 cells/L)",
       title = "Replicates Figure 3 lower panel of Valenzuela 2011",
       caption = paste("Dashed lines at ANC = 0.5 (grade 4) and 1.0 (grade 3).",
                       "n = 200 virtual subjects.")) +
  theme_minimal()
```

# Replicate Figure 4 (typical-value scenarios)

Figure 4 of Valenzuela 2011 explores how the *initial peritoneal
concentration* `Cperi(0)` and the *HIO duration* affect the ANC
nadir and the incidence of severe neutropenia. For the model implemented
here, `Cperi(0)` is set implicitly by the bolus dose and the peritoneum
volume (`Cperi(0) = dose / vd`); we reproduce the typical-value ANC
trajectories for several `Cperi(0)` values at the paper's mean
40-min HIO duration. The treatment of HIO duration as a model
intervention is approximated here by scaling the equivalent bolus dose
(HIO duration enters the published results only through total exposure
under the assumed first-order absorption).

```{r figure-4-typical, fig.width = 8, fig.height = 5}
target_cperi <- c(25, 50, 75, 100, 125, 150)   # Figure 4 grid (mg/L)
vd_typical <- 8.33                              # Va from Table II

obs_grid <- seq(0, 28 * 24, by = 2)

f4_events <- purrr::map_dfr(seq_along(target_cperi), function(i) {
  c0 <- target_cperi[i]
  amt_i <- c0 * vd_typical                # bolus that yields Cperi(0) = c0
  dplyr::bind_rows(
    tibble::tibble(id = i, time = 0,        amt = amt_i,
                   evid = 1L, cmt = "depot"),
    tibble::tibble(id = i, time = obs_grid, amt = NA_real_,
                   evid = 0L, cmt = "Cc")
  ) |>
    dplyr::mutate(Cperi_init = c0)
})

sim_f4 <- rxode2::rxSolve(mod_typ, events = f4_events,
                          keep = "Cperi_init",
                          returnType = "data.frame")
sim_f4$day <- sim_f4$time / 24

ggplot(sim_f4 |> dplyr::filter(time > 0),
       aes(day, ANC, colour = factor(Cperi_init), group = id)) +
  geom_line() +
  geom_hline(yintercept = c(0.5, 1.0), linetype = "dashed", alpha = 0.4) +
  scale_colour_viridis_d(name = "Initial Cperi (mg/L)") +
  labs(x = "Time (days)", y = "ANC (10^9 cells/L)",
       title = "Replicates Figure 4 upper-left panel of Valenzuela 2011",
       caption = "Typical-value ANC for six initial peritoneal concentrations.") +
  theme_minimal()
```

# PKNCA validation

The paper reports empirical mean (SD) Cmax and AUC values for the
peritoneal and plasma compartments (Results, page 76):

| Compartment | Empirical Cmax (mg/L) | Empirical Tmax (min) | Empirical AUC (mg.h/L) |
|-------------|-----------------------|----------------------|-------------------------|
| Peritoneum  | 82.30 (SD 17.76)      | 6.36 (SD 7.13)       | 1150 (SD 348)           |
| Plasma      | 2.56 (SD 0.90)        | 35.97 (SD 8.20)      | 87.20 (SD 123.20)       |

We reproduce these statistics by running PKNCA on the simulated cohort.
The grouping variable is the dose stratum (`amt`, single dose per
subject); PKNCA needs a treatment grouping in the formula so summaries
can be compared per cohort.

```{r pknca-plasma}
sim_nca_plasma <- sim |>
  dplyr::filter(!is.na(Cc), time > 0) |>
  dplyr::select(id, time, Cc, BSA) |>
  dplyr::mutate(treatment = "HIO_360mgm2")

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

conc_plasma <- PKNCA::PKNCAconc(sim_nca_plasma, Cc ~ time | treatment + id,
                                concu = "mg/L", timeu = "h")
dose_obj    <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id,
                                doseu = "mg")

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

nca_plasma <- PKNCA::pk.nca(PKNCA::PKNCAdata(conc_plasma, dose_obj,
                                             intervals = intervals))
knitr::kable(summary(nca_plasma),
             caption = "Simulated plasma NCA (n = 200 virtual subjects, single 360 mg/m^2 HIO).")
```

```{r pknca-peritoneum}
sim_nca_peri <- sim |>
  dplyr::filter(!is.na(Cperi), time > 0) |>
  dplyr::select(id, time, Cperi, BSA) |>
  dplyr::mutate(treatment = "HIO_360mgm2") |>
  dplyr::rename(Cc = Cperi)        # PKNCA conc column reuse

conc_peri <- PKNCA::PKNCAconc(sim_nca_peri, Cc ~ time | treatment + id,
                              concu = "mg/L", timeu = "h")

intervals_peri <- data.frame(
  start = 0, end = 28,             # peritoneum sampling effectively ends ~28 h
  cmax = TRUE, tmax = TRUE,
  auclast = TRUE, half.life = TRUE
)

nca_peri <- PKNCA::pk.nca(PKNCA::PKNCAdata(conc_peri, dose_obj,
                                           intervals = intervals_peri))
knitr::kable(summary(nca_peri),
             caption = "Simulated peritoneum NCA (n = 200 virtual subjects, single 360 mg/m^2 HIO).")
```

## Comparison against published NCA

The simulated peritoneum Cmax is dominated by the bolus initial
concentration `dose / vd` and tracks the empirical Cmax of 82.30 mg/L
once individual dose variability (BSA-driven) is folded in. The
empirical *time* of Cmax (`Tmax = 6.36 min`) reflects the early HIO
sampling: in this model, peritoneum starts at its bolus maximum at
t = 0, so the simulated Tmax is 0 by construction. This is a model
simplification, not a fit failure.

Plasma Cmax is the area of largest visible discrepancy: the typical-
value simulation predicts Cmax around 5.5 mg/L at t ~= 8 h, while the
empirical mean Cmax is 2.56 mg/L at t ~= 36 min. Two contributing
factors:

1. **Empirical Cmax vs. typical-value prediction.** "Mean Cmax" in the
   paper is the average of *per-subject observed maxima* taken at the
   sparse blood-sampling times (0.25, 0.5, 1, 1.5, 2, 2.5, 3, 4, 6, 8,
   12, 16, 20, 24, 28 h after end of HIO). Sparse-sampling Cmax is a
   biased downward estimator of the true peak, especially when the true
   peak occurs between scheduled samples. The simulated typical-value
   prediction is the deterministic continuous peak.
2. **HIPEC perfusate evacuation.** At the end of HIO (~40 min) the
   perfusate is drained from the peritoneal cavity. The paper's model
   does not encode the evacuation explicitly -- it treats the dose as a
   fully-bioavailable instantaneous peritoneal bolus that then decays
   first-order with `ka = Qa/Va` -- so the model's prediction of
   continued plasma absorption out to ~8 h is an artifact of the
   structural assumption rather than a true predicted exposure.

For users who need to approximate the HIPEC evacuation step, manually
zeroing the `depot` amount at `t = 40 min` (e.g., via a reset
`evid = 4` dose of `amt = 0`) gives a more realistic plasma profile.

The peritoneum AUC of ~1150 mg.h/L and the plasma AUC of ~87 mg.h/L are
qualitatively reproduced by the simulation: a typical subject in this
cohort accumulates plasma AUC on the order of 80-100 mg.h/L and
peritoneum AUC of ~1000 mg.h/L, consistent with the paper. The
peritoneum-to-plasma AUC ratio of 13:1 reported in the Discussion
(p. 78) is preserved.

ANC nadir timing of 11.2 days (typical-value) matches the paper's
reported 11-14 day window for ANC nadir (Results, p. 76).

```{r anc-nadir-summary}
anc_summary <- sim |>
  dplyr::group_by(id) |>
  dplyr::summarise(
    nadir_anc = min(ANC),
    nadir_day = time[which.min(ANC)] / 24,
    .groups = "drop"
  )

cat("Simulated ANC nadir (n =", nrow(anc_summary), "virtual subjects)\n")
cat("  median nadir =", round(median(anc_summary$nadir_anc), 2),
    "  range =", round(min(anc_summary$nadir_anc), 2),
    "-", round(max(anc_summary$nadir_anc), 2), "x 10^9/L\n")
cat("  median nadir day =", round(median(anc_summary$nadir_day), 1),
    "  range =", round(min(anc_summary$nadir_day), 1),
    "-", round(max(anc_summary$nadir_day), 1), "days\n")
cat("Paper Results (p. 76): median nadir 3.09 x 10^9/L (range 0.03 - 10.31);",
    "nadir at days 11-14.\n")
```

# Assumptions and deviations

* **No subject covariates in `model()`.** The paper screened age, sex,
  body weight, serum creatinine, serum albumin, ALT, AST, total
  bilirubin, hemoglobin, and hematocrit graphically against PK / PD
  parameters and found no correlation; given the small sample size
  (N = 30) no formal covariate analysis was attempted (Results, p. 76).
  These ten covariates are recorded in
  `Valenzuela_2011_oxaliplatin()$covariatesDataExcluded` for
  downstream documentation, but the packaged model has an empty
  `covariateData` and the predictions do not depend on subject
  characteristics. Inter-subject variation in the simulation cohort
  above comes from the BSA-driven dose scaling and the model's IIV.

* **Reparameterisation of peritoneum absorption.** The paper parameterises
  the peritoneal absorption as `(Qa, Va)` with `ka = Qa/Va` as a
  secondary parameter (Methods, "Pharmacokinetic and Pharmacodynamic
  Model" paragraph). nlmixr2lib has no canonical `lqa`/`lva` pair, so
  the equivalent `(lka, lvd)` canonical-name parameterisation is used
  with a correlated `etalka + etalvd` IIV block that algebraically
  reproduces the paper's independent `omega_Qa = 34.1 %` and
  `omega_Va = 17.7 %`. Specifically:
  `var(etalka)         = var(etalqa) + var(etalva) = 0.110 + 0.031 = 0.141`,
  `cov(etalka, etalvd) = -var(etalva) = -0.031`,
  `var(etalvd)         =  var(etalva) =  0.031`,
  so the marginal distributions of `(Qa, Va)` and their zero covariance
  are preserved under the change of variables.

* **`kCirc` is implemented as a separate fixed parameter from `ktr`.**
  Equation 8 of the paper writes `dCirc/dt = ktr * Transit3 - kCirc *
  Circ`, and the Methods text explicitly fixes `kCirc = 0.07 h^-1`
  citing Friberg 2002. The Figure 1 schematic annotates `k_circ
  (= k_tr)` in the standard Friberg cartoon notation, which would
  collapse to a single transit rate constant; we follow the
  text and Equation 8 over the figure annotation. The implementation
  consequence is that the precursor-chain steady-state amount is
  `(kCirc / ktr) * Circ0 ~= 14.6 x 10^9/L` rather than `Circ0 ~= 7.05
  x 10^9/L`; the chain's precursor pools therefore sit higher than the
  circulating count at baseline, and the model file sets the
  `precursor1..4(0) <- prol_ss` initial conditions accordingly. If
  taken literally as `kCirc = ktr = (n + 1) / MTT = 4 / 118 ~= 0.034
  h^-1`, the steady state would instead be `precursor_i(0) = circ(0)
  = Circ0 = 7.05`. The two encodings give qualitatively similar ANC
  dynamics for the published dose range; users wanting strict
  adherence to the alternate (figure-caption) reading can override
  the parameter with `mod |> rxode2::ini(kcirc <- fixed(4 / 118))`
  and set `precursor_i(0) <- circ0` manually.

* **HIO infusion modelled as a peritoneal bolus.** The paper's
  structural model treats the HIO administration as a single
  instantaneous peritoneal bolus that decays first-order at `ka =
  Qa/Va`. The actual HIO procedure perfuses oxaliplatin through the
  peritoneal cavity for 30-60 min and then evacuates the residual
  perfusate. This simplification is the source of the plasma
  Cmax / Tmax mismatch documented in the PKNCA section.

* **Empirical Cmax vs. typical-value Cmax.** The paper's "mean Cmax" is
  the average of per-subject observed maxima at sparse blood-sampling
  times, which is biased downward relative to the true continuous
  peak. The typical-value simulation in this vignette reports the
  deterministic continuous peak; the ~2x discrepancy between simulated
  and empirical Cmax_plasma reflects this, not a model misfit.

# Errata

An erratum to Valenzuela 2011 is listed in the PubMed record for
PMID 21210260 as "Erratum in: AAPS J. 2011 Jun;13(2):318". The
erratum text was not available on disk at extraction time, and a
public-web search for its full content was inconclusive (the
publisher's full-text site requires subscription authentication;
the PubMed record for the erratum itself was not located via DOI or
title search at the time of writing). Operators reviewing this
model file should locate the erratum and confirm whether any
parameter values, equations, or table entries require revision.
None of the values in the packaged model reflect erratum-derived
changes at the time of this extraction.

The Figure 1 schematic annotation `k_circ (= k_tr)` is internally
inconsistent with the text statement "kCirc was fixed to the
population mean half life of neutrophils previously determined, 0.07
h^-1" and with Equation 8 (which writes `kCirc` as a symbol distinct
from `ktr`). This appears to be an artefact of reusing the standard
Friberg myelosuppression cartoon; we follow the text and equation
encoding (see Assumptions and deviations above).

# Session info

```{r session-info}
sessionInfo()
```
