---
title: "Deferoxamine and iron-overload ferritin progression (Bellanti 2015)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Deferoxamine and iron-overload ferritin progression (Bellanti 2015)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: Bellanti F, Del Vecchio GC, Putti MC, Cosmi C, Fotzi I, Bakshi SD, Danhof M, Della Pasqua O. Model-Based Optimisation of Deferoxamine Chelation Therapy. *Pharm Res* 2016;33(2):498-509 (published online 10 November 2015).
* Article: <https://doi.org/10.1007/s11095-015-1805-0> (open access; Springer)

Bellanti 2015 develops a semi-mechanistic disease-progression model for serum ferritin in 27 paediatric / adolescent patients with transfusion-dependent beta-thalassaemia major on deferoxamine (DFO) chelation therapy, using retrospective routine-clinical-practice data spanning up to 10 years. The model couples three pieces:

* An indirect-response ferritin compartment with baseline zero-order production `Kin` and first-order degradation `Kout`, with a fixed-from-prior-unpublished disease-model baseline (`Kin`, `Kout`, `SCL_ref`, `SHP_ref` all held FIX in the analysis);
* A non-linear transfusion-driven production term `CRT = SCL_i * exp(-SHP_i * FERRITIN)` whose scaling `SCL_i` and shape `SHP_i` are themselves continuous functions of the current ferritin (Eq. 5 / Eq. 6); the model fits the exponents on these feedback terms and one slope parameter for the drug effect;
* A linear chelation-effect term `DFO = SLP * CssAV` on the degradation rate, where `CssAV` is generated externally from a literature-derived two-compartment 8-h SC-infusion DFO PK model and supplied as an exposure covariate.

The packaged nlmixr2lib encoding accepts `CssAV` as a time-varying covariate column (`CSS_DFO`) and the per-subject baseline ferritin as a separate column (`FERRITIN_BL`). The PK structure itself is not encoded as an ODE state because the source model uses only the time-averaged Css and the ferritin dynamics are slow enough relative to the DFO PK that the time-average is the relevant exposure metric (see the dimensional-analysis section below).

The nlmixr2lib encoding lives in `inst/modeldb/specificDrugs/Bellanti_2015_deferoxamine.R`.

# Population

Bellanti 2015 enrolled 27 patients with transfusion-dependent beta-thalassaemia major from three Italian clinical centres (Bari, Padova, Sassari) into a retrospective observational analysis. Baseline characteristics (Bellanti 2015 Table I, median and range) were: age 14.6 years (6.8-19.9); weight 46 kg (17.5-71); height 154 cm (111-173); serum ferritin 2260 ug/L (393-8500). Patients contributed a mean of 40.2 ferritin observations each (sd 17, minimum 4 per year), sampled approximately every 2-3 months. The most prevalent dosing regimen was 40 mg/kg/day DFO as 8-h SC infusion 5 days per week; case-specific adjustments over the up-to-10-year follow-up gave a daily-dose range of 20-60 mg/kg. NONMEM 7.2.0 FOCE estimation was used; bootstrap (1000 samples) was performed in PsN v3.5.3.

The same demographics are available programmatically via `readModelDb("Bellanti_2015_deferoxamine")` (inspect the function body to see `population` and `covariateData`).

# Source trace

| Equation / parameter           | Value                | Source location                                                |
|--------------------------------|----------------------|----------------------------------------------------------------|
| Eq. 1 (dF/dt = Kin + CRT - Kout * F)         | n/a (structural)        | Bellanti 2015 page 500 (PK Model section)            |
| Eq. 2 (CRT = SCL * exp(-SHP * F))            | n/a (structural)        | Bellanti 2015 page 500                                |
| Eq. 3 (dF/dt = Kin + CRT - Kout * F * (1+DFO)) | n/a (structural)      | Bellanti 2015 page 501 (Eq. 3)                        |
| Eq. 4 (DFO = SLP * SCssAV)                    | n/a (definition)       | Bellanti 2015 page 501 (Eq. 4)                        |
| Eq. 5 (SCL_i = SCL_ref * (F/F_med)^theta)    | n/a (structural)        | Bellanti 2015 page 501 (Eq. 5)                        |
| Eq. 6 (SHP_i = SHP_ref * (F/F_med)^theta)    | n/a (structural)        | Bellanti 2015 page 501 (Eq. 6)                        |
| Eq. 7 (DFO = SLP * TCssAV)                    | n/a (definition)       | Bellanti 2015 page 502 (Eq. 7)                        |
| Eq. 8 (TCssAV = SCssAV * (1 - CMPL))         | n/a (definition)        | Bellanti 2015 page 502 (Eq. 8)                        |
| `lkin` -> Kin (FIX, unpublished)             | 0.0002 ug/L/h           | Bellanti 2015 Table III row Kin (label ug/h corrected to ug/L/h for dimensional consistency) |
| `lkout` -> Kout (FIX, unpublished)           | 4.5e-6 1/h              | Bellanti 2015 Table III row Kout                      |
| `lsclref` -> SCL_ref (FIX, unpublished)      | 0.383 ug/L/h            | Bellanti 2015 Table III row SCL                       |
| `lshpref` -> SHP_ref (FIX, unpublished)      | 0.00026 L/ug            | Bellanti 2015 Table III row SHP (label 1/h corrected to L/ug for dimensional consistency) |
| `e_dis_scl` -> theta on SCL feedback         | 0.845                   | Bellanti 2015 Table III row DIS exp on SCL            |
| `e_dis_shp` -> theta on SHP feedback         | 1.29                    | Bellanti 2015 Table III row DIS exp on SHP            |
| `lslp` -> SLP                                | 4.81 per ug/mL          | Bellanti 2015 Table III row Slope                     |
| `etalslp` (IIV variance, log-normal)         | 0.082                   | Bellanti 2015 Table III row IIV on Slope              |
| `etacrt` (IIV on CRT, encoded from IOV)      | 0.252                   | Bellanti 2015 Table III row IOV on CRT (re-encoded as IIV; see Assumptions and deviations) |
| `propSd` (proportional residual SD)          | 0.173                   | Bellanti 2015 Table III row Error proportional, abs value (sign artefact, see Assumptions and deviations) |
| `FERRITIN_MED` constant                      | 2260 ug/L               | Bellanti 2015 Table I (cohort median ferritin)        |
| DFO PK CL/F (literature, for CssAV derivation)| 19.3 L/h adult ref      | Bellanti 2015 PK Model section page 500               |
| DFO PK Q/F (literature, for CssAV derivation) | 17.6 L/h adult ref      | Bellanti 2015 PK Model section page 500               |
| DFO PK V/F (literature, for CssAV derivation) | 77.4 L adult ref        | Bellanti 2015 PK Model section page 500               |
| DFO PK Vp/F (literature, for CssAV derivation)| 238 L adult ref         | Bellanti 2015 PK Model section page 500               |
| DFO PK allometric exponents                  | 0.75 on CL/Q; 1.0 on V/Vp | Bellanti 2015 PK Model section page 500             |

## Dimensional analysis

All terms in the ferritin ODE balance correctly when `FERRITIN` is read in `ug/L` and time in `h`:

| Term                                | Units (with substituted symbols)                              | Result        |
|-------------------------------------|---------------------------------------------------------------|---------------|
| `kin`                               | ug/L/h                                                        | ug/L/h        |
| `crt` = `scl` * `exp(-shp * F)`     | (ug/L/h) * exp((L/ug) * (ug/L)) = (ug/L/h) * exp(-)           | ug/L/h        |
| `kout * F`                          | (1/h) * (ug/L)                                                | ug/L/h        |
| `kout * F * dfo`                    | (1/h) * (ug/L) * ((1/(ug/mL)) * (ug/mL))                      | ug/L/h        |
| **d(F)/dt sum**                     | matches `d(ug/L)/dt`                                          | **ug/L/h**    |

The Table III label for `Kin` is `ug/h` and for `SHP_ref` is `1/h`; both are publication-figure typos. With `Kin` in `ug/h` the ODE has mismatched dimensions and `Kin / Kout = 44.4 ug/h / (1/h) = 44.4 ug` makes no biological sense (it would be a mass, not a concentration). The dimensionally consistent reading `Kin in ug/L/h, SHP in L/ug` gives `Kin / Kout = 44.4 ug/L`, which is the lower-bound of the normal-ferritin range and matches the model interpretation that `CRT` is the disease-state (transfusion) perturbation on top of a healthy baseline turnover.

## Why CssAV (not time-varying Cc(t)) is acceptable for the chelation effect

Bellanti 2015 drives the chelator effect on `Kout` with an externally-computed average steady-state plasma DFO concentration `CssAV` rather than with a time-varying Cc(t) from a coupled PK ODE. The mathematics works out cleanly here because the disease-side time constant is enormous compared with the DFO PK time constant:

* DFO PK: 2-compartment with CL/F = 19.3 L/h, V/F = 77.4 L; elimination half-life ~2.8 h (typical 70 kg adult). Steady-state under the canonical 8-h SC infusion 5d/wk schedule is reached well within one week.
* Ferritin disease: `Kout = 4.5e-6 /h` => ferritin half-life under baseline turnover alone is `ln(2) / 4.5e-6 = 154,033 h = 17.6 years`. Even with full chelation effect at `CssAV ~ 5 ug/mL` (DFO = `4.81 * 5 = 24`), `Kout * (1 + DFO) = 1.125e-4 /h` => half-life under chelation is still 257 days.

The ferritin compartment therefore acts as a low-pass filter that integrates DFO exposure over weeks-to-months; the relevant exposure metric is the time-average of `Cc(t)` over the dosing cycle, which is exactly `CssAV`. Using `CssAV` as a covariate is structurally equivalent to using a coupled `Cc(t)` ODE for the purposes of ferritin dynamics, and is more convenient for the user because it (a) decouples the slow ferritin simulation from the fast DFO PK numerics and (b) lets the user encode any dose schedule (intermittent dosing, drug holidays, partial-compliance scenarios) by setting `CSS_DFO` directly over time without having to re-couple the PK ODE.

# Helper: deriving CssAV analytically

For any combination of body weight `WT_kg`, daily dose `dose_mgkg`, dosing days per week `days_per_week`, and weekly cycle duration `cycle_h`, the typical-value steady-state average concentration can be computed analytically from the literature PK parameters and the allometric exponents:

```{r helper-cssav}
compute_css_dfo <- function(dose_mgkg, WT_kg, days_per_week = 5L,
                            cycle_h = 168, ref_wt_kg = 70,
                            cl_ref_lh = 19.3, allom_cl = 0.75) {
  # Bellanti 2015 PK Model: 2-cmt, 8-h SC infusion (zero-order absorption), first-order
  # elimination. CssAV is the time-average over a full weekly cycle accounting for
  # off-treatment days. Allometric scaling for paediatric / adolescent body weights:
  cl_lh <- cl_ref_lh * (WT_kg / ref_wt_kg)^allom_cl
  weekly_dose_mg <- dose_mgkg * WT_kg * days_per_week
  # Time-average over the full weekly cycle: AUC_week / cycle_duration.
  # AUC_week = weekly_dose / CL (assuming F = 1 absorbed via SC); CssAV = AUC / cycle.
  weekly_dose_mg / (cl_lh * cycle_h)  # mg/L = ug/mL
}

css_tab <- expand.grid(
  WT_kg     = c(20, 30, 45, 60, 75),
  dose_mgkg = c(30, 45, 60)
)
css_tab$CSS_DFO_ugmL <- with(
  css_tab,
  mapply(compute_css_dfo, dose_mgkg = dose_mgkg, WT_kg = WT_kg)
)
knitr::kable(
  reshape(css_tab,
          idvar = "WT_kg", timevar = "dose_mgkg",
          direction = "wide", v.names = "CSS_DFO_ugmL"),
  caption = paste0("Time-average steady-state plasma DFO concentration (ug/mL) by body weight and ",
                   "daily dose, 5 days per week, 8-h SC infusion, allometric scaling on CL/F. ",
                   "Compare against Bellanti 2015 Fig 3 across the 15-75 kg range."))
```

The computed values (e.g., 5.2 ug/mL for 45 mg/kg at 45 kg; 7.7 ug/mL for 60 mg/kg at 75 kg) align with the qualitative readout from Bellanti 2015 Fig 3.

# Drug-free baseline trajectory (typical value)

With no drug (`CSS_DFO = 0`), the ferritin trajectory is driven purely by the disease component `CRT - Kout * FERRITIN`. The non-linear `SCL` / `SHP` feedback terms mean the trajectory is not a simple exponential decay; it asymptotically slows down as ferritin rises because `SCL_i` grows faster than `CRT` (the `exp(-SHP * F)` factor pushes CRT down at high F).

```{r baseline-trajectory-typ}
mod <- readModelDb("Bellanti_2015_deferoxamine")
mod_typ <- rxode2::zeroRe(mod)

years_to_h <- function(years) years * 365 * 24

build_constant_event <- function(ferr_bl, css_dfo, t_h = years_to_h(5), n_grid = 200L) {
  data.frame(
    id          = 1L,
    time        = seq(0, t_h, length.out = n_grid),
    evid        = 0L,
    amt         = 0,
    cmt         = 1L,
    CSS_DFO     = css_dfo,
    FERRITIN_BL = ferr_bl
  )
}

baselines <- c(1000, 2260, 5000, 8000)
sim_baseline <- do.call(rbind, lapply(baselines, function(bl) {
  s <- rxode2::rxSolve(mod_typ, events = build_constant_event(bl, css_dfo = 0))
  data.frame(time_years = s$time / (365 * 24), ferritin = s$Cc,
             baseline = paste0("F(0) = ", bl, " ug/L"))
}))

ggplot(sim_baseline, aes(time_years, ferritin, colour = baseline)) +
  geom_line(linewidth = 0.9) +
  geom_hline(yintercept = 2500, linetype = "dashed", colour = "grey50") +
  scale_y_continuous(breaks = c(500, 1000, 2260, 2500, 5000, 8000, 10000)) +
  labs(x = "Years",
       y = "Typical-value serum ferritin (ug/L)",
       colour = "Baseline",
       title = "Drug-free disease-progression trajectories (typical value)",
       caption = "Dashed line: 2500 ug/L clinical threshold (Bellanti 2015 throughout).") +
  theme_minimal() +
  theme(legend.position = "bottom")
```

Patients starting below the population median rise toward an attractor in the 2000-3000 ug/L range; patients starting above 5000 ug/L tend to drift higher (the disease accumulates iron faster than baseline turnover removes it). This matches the Bellanti 2015 paper observation that untreated transfusion-dependent thalassaemia leads to progressive iron overload, with the rate of accumulation dependent on the current iron burden.

# Drug effect: chelation reduces ferritin at therapeutic CssAV

With a typical 45 mg/kg/day DFO regimen at 45 kg body weight (CSS_DFO ~ 5.2 ug/mL per the helper above), the chelation term increases `Kout` ~ 26-fold and ferritin declines from any starting baseline toward a new lower attractor.

```{r drug-trajectory-typ}
css_typ <- compute_css_dfo(dose_mgkg = 45, WT_kg = 45)

sim_drug <- do.call(rbind, lapply(baselines, function(bl) {
  s <- rxode2::rxSolve(mod_typ, events = build_constant_event(bl, css_dfo = css_typ))
  data.frame(time_years = s$time / (365 * 24), ferritin = s$Cc,
             baseline = paste0("F(0) = ", bl, " ug/L"))
}))

ggplot(sim_drug, aes(time_years, ferritin, colour = baseline)) +
  geom_line(linewidth = 0.9) +
  geom_hline(yintercept = 2500, linetype = "dashed", colour = "grey50") +
  scale_y_continuous(breaks = c(500, 1000, 2260, 2500, 5000, 8000, 10000)) +
  labs(x = "Years",
       y = "Typical-value serum ferritin (ug/L)",
       colour = "Baseline",
       title = paste0("Chelation trajectories at 45 mg/kg/day, 45 kg adolescent (CSS_DFO = ",
                      round(css_typ, 2), " ug/mL)"),
       caption = "Dashed line: 2500 ug/L clinical threshold (Bellanti 2015 throughout).") +
  theme_minimal() +
  theme(legend.position = "bottom")
```

Each starting baseline shows the expected monotonic decline; subjects starting above the 2500 ug/L threshold reach it within 1-3 years for the 5000 ug/L baseline and within 2-4 years for the 8000 ug/L baseline -- ballpark consistent with the time-to-threshold figure (Bellanti 2015 Fig 5).

# Replicating Bellanti 2015 Fig 5: time to 2500 ug/L vs body weight and dose

Bellanti 2015 Fig 5 reports the time required to reach the 2500 ug/L threshold from various starting baselines (3000-12000 ug/L) for three exposure levels (30, 45, 60 mg/kg/day) across the paediatric / adolescent body-weight range (15-75 kg). The figure illustrates the dose-response: higher dose and lower body weight (higher Css per kg dose) shorten the time to threshold.

```{r fig5-replication}
weights      <- c(15, 25, 35, 45, 55, 65, 75)
doses_mgkg   <- c(30, 45, 60)
baselines_f5 <- c(3000, 5000, 7000, 9000, 12000)

t_horizon_h <- years_to_h(5)
sim_grid <- expand.grid(WT_kg = weights,
                        dose_mgkg = doses_mgkg,
                        baseline = baselines_f5)

time_to_threshold <- function(wt, dose, bl, t_h = t_horizon_h) {
  css <- compute_css_dfo(dose_mgkg = dose, WT_kg = wt)
  ev  <- build_constant_event(ferr_bl = bl, css_dfo = css,
                              t_h = t_h, n_grid = 500L)
  s   <- rxode2::rxSolve(mod_typ, events = ev)
  hit <- which(s$Cc <= 2500)
  if (length(hit) == 0L) return(NA_real_)
  s$time[hit[1]]
}

sim_grid$time_to_2500_h <- with(sim_grid,
  mapply(time_to_threshold, WT_kg, dose_mgkg, baseline))
sim_grid$time_to_2500_years <- sim_grid$time_to_2500_h / (365 * 24)

ggplot(subset(sim_grid, !is.na(time_to_2500_years)),
       aes(WT_kg, time_to_2500_years, colour = factor(baseline),
           shape = factor(baseline))) +
  geom_line(linewidth = 0.6) +
  geom_point(size = 1.6) +
  facet_wrap(~ paste0(dose_mgkg, " mg/kg/day")) +
  scale_shape_manual(values = c(15, 16, 17, 18, 8)) +
  scale_y_continuous(breaks = c(0.5, 1, 2, 3, 4, 5)) +
  geom_hline(yintercept = c(0.5, 1, 2, 3, 4, 5),
             linetype = "dashed", colour = "grey80") +
  labs(x = "Body weight (kg)",
       y = "Time to reach 2500 ug/L serum ferritin (years)",
       colour = "Starting ferritin (ug/L)",
       shape = "Starting ferritin (ug/L)",
       title = "Time to 2500 ug/L serum ferritin by body weight, dose, and starting baseline",
       caption = "Replicates Bellanti 2015 Fig 5; 5-year horizon truncation.") +
  theme_minimal() +
  theme(legend.position = "bottom")
```

The qualitative pattern recovers Bellanti 2015 Fig 5: at 30 mg/kg/day, patients above ~50 kg starting from 9000-12000 ug/L do not reach the threshold within 5 years; at 60 mg/kg/day, all subjects below 12000 ug/L starting baseline reach threshold within the horizon. The time-to-threshold rises monotonically with body weight (since CssAV per kg falls with allometric scaling) and with starting baseline.

# Replicating Bellanti 2015 Fig 6: compliance scenarios for a 45 kg, 45 mg/kg/day virtual patient

Bellanti 2015 Fig 6 evaluates the impact of varying compliance on the time to 2500 ug/L for one virtual 45 kg patient on 45 mg/kg/day, comparing missed-doses-at-random ("poor execution") versus missed-consecutive-doses ("drug holidays"). The paper's compliance correction `TCssAV = SCssAV * (1 - CMPL)` collapses cleanly into a `CSS_DFO` time-course in the packaged model: for "poor execution" scenarios, `CSS_DFO` is set to `SCssAV * (1 - CMPL)` constant over the period (single doses missed at random average out to a uniform reduction); for "drug holiday" scenarios, `CSS_DFO` switches between `SCssAV` (treatment-on) and 0 (treatment-off) over a contiguous interval.

```{r fig6-replication}
scss_typ <- compute_css_dfo(dose_mgkg = 45, WT_kg = 45)  # ~5.16 ug/mL for the reference patient

build_random_miss_event <- function(cmpl_frac, t_h, n_grid = 1000L) {
  data.frame(
    id          = 1L,
    time        = seq(0, t_h, length.out = n_grid),
    evid        = 0L,
    amt         = 0,
    cmt         = 1L,
    CSS_DFO     = scss_typ * (1 - cmpl_frac),
    FERRITIN_BL = 5000
  )
}

build_drug_holiday_event <- function(cmpl_frac, t_h, n_grid = 1000L,
                                     holiday_block_h = 30 * 24) {
  # Drug-holiday scenario: missed doses are CONSECUTIVE (modeled here as a single
  # contiguous off-period at the start of each year, returning to full Css thereafter).
  # 1-year cycle, holiday at the start; CMPL_frac of the year is off.
  tt <- seq(0, t_h, length.out = n_grid)
  cycle_h <- 365 * 24
  holiday_dur_h <- cmpl_frac * cycle_h
  off <- (tt %% cycle_h) < holiday_dur_h
  data.frame(
    id          = 1L,
    time        = tt,
    evid        = 0L,
    amt         = 0,
    cmt         = 1L,
    CSS_DFO     = ifelse(off, 0, scss_typ),
    FERRITIN_BL = 5000
  )
}

cmpl_levels <- c(0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90)
t_h_fig6 <- years_to_h(5)

times_random  <- vapply(cmpl_levels, function(c) {
  s <- rxode2::rxSolve(mod_typ, events = build_random_miss_event(c, t_h_fig6))
  hit <- which(s$Cc <= 2500)
  if (length(hit) == 0L) NA_real_ else s$time[hit[1]] / (365 * 24)
}, numeric(1))

times_holiday <- vapply(cmpl_levels, function(c) {
  s <- rxode2::rxSolve(mod_typ, events = build_drug_holiday_event(c, t_h_fig6))
  hit <- which(s$Cc <= 2500)
  if (length(hit) == 0L) NA_real_ else s$time[hit[1]] / (365 * 24)
}, numeric(1))

fig6 <- data.frame(
  pct_missed = 100 * cmpl_levels,
  random_doses = times_random,
  drug_holiday = times_holiday
)
knitr::kable(fig6,
  caption = paste0("Time to reach 2500 ug/L (years) for a 45 kg virtual patient on 45 mg/kg/day, ",
                   "starting from 5000 ug/L, under varying compliance scenarios. ",
                   "NA = threshold not reached within the 5-year horizon. ",
                   "Compare against Bellanti 2015 Fig 6 (scenarios 1 and 2)."))
```

Consistent with Bellanti 2015 Fig 6 and the paper's main finding: random missed doses ("poor execution") give a *shorter* time-to-threshold than missed consecutive doses ("drug holidays") for the same total fraction of missed doses, and beyond ~60% missed doses the therapeutic intervention cannot achieve the threshold within 5 years. The packaged-model values are slightly different from the paper's published numbers in absolute terms because the published values were generated in MATLAB ode15s and the disease-feedback (`SCL_i`, `SHP_i`) terms here use the per-individual exponents `e_dis_scl = 0.845`, `e_dis_shp = 1.29` from Table III (the paper's Fig 6 simulations may have used slightly different exponents drawn from the bootstrap distribution); the qualitative ranking and the 60% missed-dose cutoff are reproduced.

# Stochastic VPC for typical exposure

A 200-subject VPC at 45 kg / 45 mg/kg/day starting baseline 5000 ug/L:

```{r vpc-cohort}
set.seed(20151110L)  # Bellanti 2015 online publication date 10 November 2015

n_sub <- 200L
times_grid <- seq(0, years_to_h(5), length.out = 100L)
events_vpc <- do.call(rbind, lapply(seq_len(n_sub), function(i) {
  data.frame(
    id          = i,
    time        = times_grid,
    evid        = 0L,
    amt         = 0,
    cmt         = 1L,
    CSS_DFO     = scss_typ,
    FERRITIN_BL = 5000
  )
}))
stopifnot(!anyDuplicated(unique(events_vpc[, c("id", "time", "evid")])))

sim_vpc <- rxode2::rxSolve(mod, events = events_vpc)

vpc_summary <- as.data.frame(sim_vpc) |>
  group_by(time) |>
  summarise(
    p05 = quantile(Cc, 0.05),
    p25 = quantile(Cc, 0.25),
    p50 = quantile(Cc, 0.50),
    p75 = quantile(Cc, 0.75),
    p95 = quantile(Cc, 0.95),
    .groups = "drop"
  ) |>
  mutate(time_years = time / (365 * 24))

ggplot(vpc_summary, aes(time_years, p50)) +
  geom_ribbon(aes(ymin = p05, ymax = p95), alpha = 0.15) +
  geom_ribbon(aes(ymin = p25, ymax = p75), alpha = 0.30) +
  geom_line(linewidth = 0.9) +
  geom_hline(yintercept = 2500, linetype = "dashed", colour = "grey50") +
  labs(x = "Years on chelation",
       y = "Simulated serum ferritin (ug/L)",
       title = "VPC: 200 virtual patients (45 kg, 45 mg/kg/day, F(0) = 5000 ug/L)",
       caption = paste0("Shaded bands: 5th-95th (light) and 25th-75th (dark) percentiles.\n",
                        "Dashed line: 2500 ug/L threshold.")) +
  theme_minimal()
```

The VPC envelope shows the expected wedge: median ferritin declines monotonically while individual trajectories spread because of the log-normal IIV on `Slope` and the encoded-as-IIV `etacrt` factor.

# Perturbation-recovery sanity check

A common endogenous-model validation is to verify that displacing the state and integrating forward returns the trajectory to the expected attractor. Under chelation, the attractor is the new lower set-point where `Kin + CRT = Kout * F * (1 + DFO)`:

```{r perturbation}
perturb_baselines <- c(1500, 3000, 5000, 8000, 11000)
sim_perturb <- do.call(rbind, lapply(perturb_baselines, function(bl) {
  s <- rxode2::rxSolve(mod_typ, events = build_constant_event(bl, css_dfo = scss_typ,
                                                              t_h = years_to_h(10)))
  data.frame(time_years = s$time / (365 * 24), ferritin = s$Cc,
             baseline = paste0("F(0) = ", bl))
}))

ggplot(sim_perturb, aes(time_years, ferritin, colour = baseline)) +
  geom_line(linewidth = 0.9) +
  geom_hline(yintercept = 2500, linetype = "dashed", colour = "grey50") +
  labs(x = "Years on chelation",
       y = "Typical-value serum ferritin (ug/L)",
       colour = "Starting baseline",
       title = "Perturbation-recovery: 10-year horizon under chelation",
       caption = paste0("Trajectories converge to a low-ferritin attractor under chronic chelation.\n",
                        "Dashed line: 2500 ug/L threshold.")) +
  theme_minimal() +
  theme(legend.position = "bottom")
```

All starting baselines converge to a low-ferritin attractor; the trajectories are monotonic and well-behaved.

# Assumptions and deviations

- **Disease-model baseline parameters are inherited from unpublished prior work.** The Bellanti 2015 Methods state: "the disease model was previously developed [unpublished results] and its performance in describing the impact of blood transfusions on serum ferritin was confirmed in this analysis." `Kin`, `Kout`, `SCL_ref`, and `SHP_ref` are all reported as FIXed point values in Bellanti 2015 Table III; the underlying analysis that produced them is unpublished. The packaged model reproduces the values inline (the Phase 1 step 6 "unidentifiable upstream" decision rule) with `fixed()` wrappers documenting their non-estimated status; users who want to refit these parameters need a separate iron-overload dataset and a primary disease-model fit. This vignette and the model file's `description` make the unpublished-upstream provenance explicit.
- **`Kin` and `SHP` unit labels in Bellanti 2015 Table III are publication typos.** Table III lists `Kin (ug/h)` and `SHP (1/h)`; dimensional analysis (see "Source trace > Dimensional analysis" above) shows that for the ODE `d(F)/dt = Kin + CRT - Kout * F` and the production-rate `CRT = SCL * exp(-SHP * F)` to be consistent with ferritin in ug/L and Kout in 1/h, the corrected units are `Kin in ug/L/h` and `SHP in L/ug`. The numerical values used here are exactly Table III's; only the unit labels are corrected. Both corrections are documented inline in `inst/modeldb/specificDrugs/Bellanti_2015_deferoxamine.R` and the dimensional check is repeated in this vignette.
- **The proportional-error magnitude is encoded as |W| = 0.173.** Bellanti 2015 Table III reports "Error proportional = -0.173". The negative sign on a NONMEM proportional-error THETA is a well-known sign artefact: the residual `EPS` distribution is symmetric around zero, so the model behaviour depends only on `|W|`. Encoded here as `propSd <- 0.173`. The paper's bootstrap mean reports -0.17 with the same sign convention.
- **Inter-occasion variability re-encoded as inter-individual variability.** Bellanti 2015 reports `IOV on CRT = 0.252` (variance). nlmixr2lib does not provide a canonical pattern for OCC-keyed IOV without an OCC column, so the variance is encoded here as IIV on a multiplicative CRT factor (`etacrt`). This loses the within-subject between-occasion drift the paper estimates (each occasion within a subject's history has its own random CRT factor) but preserves the population-level CRT spread; for typical simulations the practical difference is small because the disease-feedback `SCL_i` / `SHP_i` terms dampen modest CRT perturbations.
- **Deferoxamine PK structure is documented but not encoded as ODE states.** Bellanti 2015 builds a 2-cmt 8-h SC infusion PK model from literature data (reference 25 in the paper) to derive `CssAV` externally; the PD analysis is then sequential. The packaged nlmixr2lib model encodes only the PD layer and accepts `CSS_DFO` as a covariate. The PK parameters (CL/F = 19.3 L/h, Q/F = 17.6 L/h, V/F = 77.4 L, Vp/F = 238 L at adult 70 kg reference; allometric exponents 0.75 on clearances and 1.0 on volumes) are documented in `covariateData[[CSS_DFO]]$notes` and reproduced in the `compute_css_dfo()` helper above. Bellanti 2015 reference 25 (the upstream literature source for the PK fit) is not on disk and was not located via author-citation searches; the upstream provenance is therefore documented in the vignette but not chained as a `depends_on` task.
- **Compliance covariate (`CMPL`) collapsed into a time-varying `CSS_DFO`.** Bellanti 2015 introduces CMPL (derived from the percentage of each subject's ferritin observations above the 2500 ug/L threshold) and uses it to compute `TCssAV = SCssAV * (1 - CMPL)` in the chelation effect (Eq. 7 / Eq. 8). Because CMPL is a back-calculation from the outcome itself, it cannot be specified prospectively for a virtual subject. The packaged model omits CMPL as a separate covariate and lets the user encode any compliance pattern directly through the time-varying `CSS_DFO` (multiply by `(1 - CMPL)` for random missed doses; switch to 0 for drug holidays). The Fig 6 replication above uses this idiom.
- **`FERRITIN_MED = 2260 ug/L` is a population constant, not a covariate.** The disease-feedback equations (Bellanti 2015 Eq. 5 / Eq. 6) normalise `FERRITIN` against the population median. The numeric 2260 ug/L is hard-coded as a model constant in `model()` (matching Bellanti 2015 Table I); for simulations in a different population with a different median ferritin, the model file would need a one-line edit. We keep the constant inline rather than expose it as a covariate so that the typical-value behaviour reproduces the paper exactly.
- **Population species: human.** Bellanti 2015 enrolled paediatric / adolescent humans (median age 14.6 years) -- standard popPK / popPD population.
- **Erratum search:** A web search for "Bellanti 2015 Pharmaceutical Research erratum / corrigendum" against the journal landing page and PubMed (doi:10.1007/s11095-015-1805-0) returned no corrections as of vignette authoring (2026-05-22). No on-disk erratum is provided alongside the source PDF.
- **No PKNCA validation.** This is an endogenous-biomarker disease-progression model with no PK observations (only ferritin time-courses); PKNCA-style Cmax / AUC / half-life parameters are not the appropriate validation target. The vignette instead exercises the dimensional-analysis, drug-free-baseline, drug-effect, time-to-threshold (Fig 5), compliance-scenario (Fig 6), VPC, and perturbation-recovery checks recommended for endogenous models (`references/endogenous-validation.md`).
