Skip to contents

Deferoxamine and iron-overload ferritin progression (Bellanti 2015)

Source: vignettes/articles/Bellanti_2015_deferoxamine.Rmd
Bellanti_2015_deferoxamine.Rmd

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:

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."))
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.
WT_kg CSS_DFO_ugmL.30 CSS_DFO_ugmL.45 CSS_DFO_ugmL.60
20 2.367586 3.551379 4.735171
30 2.620164 3.930246 5.240328
45 2.899688 4.349533 5.799377
60 3.115918 4.673877 6.231836
75 3.294682 4.942023 6.589364

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).

mod <- readModelDb("Bellanti_2015_deferoxamine")
mod_typ <- rxode2::zeroRe(mod)
#> ℹ parameter labels from comments will be replaced by 'label()'

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"))
}))
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'

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.

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"))
}))
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'

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.

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))
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
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.

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))
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'

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))
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'

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)."))
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).
pct_missed random_doses drug_holiday
10 1.971972 1.806807
20 NA 2.722723
30 NA 3.948949
40 NA NA
50 NA NA
60 NA NA
70 NA NA
80 NA NA
90 NA NA

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:

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)
#> ℹ parameter labels from comments will be replaced by 'label()'

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):

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))
}))
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'
#> ℹ omega/sigma items treated as zero: 'etalslp', 'etalcrt_mult'

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).