---
title: "Liposomal daunorubicin (Hempel 2003)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Liposomal daunorubicin (Hempel 2003)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

```{r meta}
mod_meta <- nlmixr2est::nlmixr(readModelDb("Hempel_2003_daunorubicin_liposomal"))$meta
```

* Citation: `r mod_meta$reference`
* Description: `r mod_meta$description`
* Article (DOI): <https://doi.org/10.1046/j.1365-2125.2003.01886.x>

This vignette validates the packaged `Hempel_2003_daunorubicin_liposomal`
model -- a one-compartment IV-infusion population PK model for total
daunorubicin (free plus liposome-encapsulated) following liposomal
daunorubicin (Daunoxome) in 24 paediatric and adolescent oncology
patients -- against the source publication's Table 2 (population PK model
development), Table 3 (BSA-normalised comparison with the prior
literature), and the representative subject profile shown in Figure 5.

# Population

Hempel 2003 pooled 214 plasma daunorubicin concentrations from 72
treatment cycles (mean 3 samples per cycle, 9 per patient) in 24
paediatric and adolescent patients enrolled across 16 German paediatric
haematology / oncology centres in the AML-REZ 97 protocol of the German
Society for Paediatric Oncology and Haematology (GPOH). Nineteen of 24
patients had relapsed acute myeloic leukaemia (AML); the remaining five
had other relapsed malignancies (two acute lymphoblastic leukaemia, one
osteosarcoma, one Ewing sarcoma, one multiple endocrine neoplasia). All
patients had previously received conventional daunorubicin and / or
doxorubicin. Per Table 1, median age was 15.4 years (range 2.84-23.2);
median body weight 48.8 kg (range 14-76.5); median height 1.61 m (range
0.89-1.98); median body surface area 1.48 m^2 (range 0.58-1.98).
Daunoxome was administered as a 1- to 2.5-hour IV infusion at 30 or 60
mg/m^2 on days 1 and 5 of each induction cycle (a pilot phase at 30
mg/m^2 was followed by escalation to 60 mg/m^2; four patients on the 60
mg/m^2 schedule intensified to days 1, 3, and 5). Total daunorubicin
(free plus liposome-encapsulated) was quantified by capillary
electrophoresis with laser-induced fluorescence detection following
organic-solvent sample preparation that destroyed the liposomes;
within-day accuracy / precision were 3.9-12.9% and 3.8-10.7% (n = 8),
and the lower limit of quantification was 2 ug/L (= 2 ng/mL).

The same information is available programmatically via the model's
`population` metadata:

```{r population}
str(mod_meta$population)
```

# Source trace

The per-parameter origin is recorded as an in-file comment next to each
`ini()` entry in
`inst/modeldb/specificDrugs/Hempel_2003_daunorubicin_liposomal.R`. The
table below collects them in one place; values come from Hempel 2003
Table 2 model 15 (the final model) unless otherwise noted.

| Parameter / equation | Value | Source location |
|----------------------|-------|-----------------|
| `lcl` (Clearance per kg body weight) | `log(0.00641)` L/h/kg | Table 2 model 15 CL = 6.41 mL/h/kg; abstract identical |
| `lvc` (Volume of distribution per kg body weight) | `log(0.0654)` L/kg | Table 2 model 15 V = 65.4 mL/kg; abstract identical |
| `etalcl` (IIV on CL) | 0.23133 | Table 2 model 15 omega_CL 51% CV; log(1 + 0.51^2) |
| `etalvc` (IIV on V)  | 0.07033 | Table 2 model 15 omega_V  27% CV; log(1 + 0.27^2) |
| `propSd` (proportional residual error) | 0.22 | Table 2 model 15 residual error 22% |
| `cl <- exp(lcl + etalcl) * WT` | n/a | Table 2 footnote *** (CL = q1 * Weight; V = q2 * Weight) |
| `vc <- exp(lvc + etalvc) * WT` | n/a | Table 2 footnote *** (CL = q1 * Weight; V = q2 * Weight) |
| `d/dt(central) <- -kel * central` | n/a | Methods (one-compartment ADVAN1 TRAN2); Results paragraph "A one-compartment model described the data adequately" |
| `Cc <- central / vc * 1000` | n/a | Unit conversion mg/L -> ng/mL (= ug/L, the paper's reported concentration unit) |
| `Cc ~ prop(propSd)` | n/a | Results paragraph "For the residual error, a proportional error model was chosen" |

Inter-occasion variability on CL of 16.7% CV (Table 2 model 15) is
explicitly noted in the source paper but is NOT encoded in the packaged
model; see Assumptions and deviations.

# Virtual cohort

The original observed daunorubicin concentrations are not publicly
available. The virtual cohort below approximates the Hempel 2003 Table 1
demographics: 24 paediatric and adolescent subjects with body weight,
height, and body surface area drawn to span the published ranges and
medians. Each subject receives a single 60 mg/m^2 IV infusion of
Daunoxome over 1.5 h (the cycle's day-1 dose); sampling follows the
schedule suggested in Methods (end of infusion, 2-4 h, 4-8 h, 12-17 h,
20-28 h post-infusion start) plus a 48-h late sample so PKNCA can
estimate the terminal-phase parameters.

```{r cohort}
set.seed(20260604L)

n_subj <- 24L

# Body weight: log-normal centred on the reported median 48.8 kg with
# the SD chosen so the simulated range covers Table 1's 14-76.5 kg.
wt_draw <- function(n) {
  s <- exp(rnorm(n, mean = log(48.8), sd = log(76.5 / 14) / 4))
  pmin(pmax(s, 14), 76.5)
}

# Body height: log-normal centred on the reported median 1.61 m with
# the SD chosen so the simulated range covers Table 1's 0.89-1.98 m.
ht_draw <- function(n) {
  s <- exp(rnorm(n, mean = log(1.61), sd = log(1.98 / 0.89) / 4))
  pmin(pmax(s, 0.89), 1.98)
}

# Body surface area via Mosteller: BSA(m^2) = sqrt(WT(kg) * HT(cm) / 3600).
# HT here is in metres, so multiply by 100 inside the square root.
bsa_mosteller <- function(wt_kg, ht_m) {
  sqrt(wt_kg * ht_m * 100 / 3600)
}

cov_tab <- tibble::tibble(
  id   = seq_len(n_subj),
  WT   = wt_draw(n_subj),
  HT   = ht_draw(n_subj)
) |>
  mutate(BSA = bsa_mosteller(WT, HT))

# Sanity check: simulated medians close to Table 1 (within ~5%).
stopifnot(
  abs(median(cov_tab$WT)  / 48.8 - 1) < 0.15,
  abs(median(cov_tab$HT)  / 1.61 - 1) < 0.15,
  abs(median(cov_tab$BSA) / 1.48 - 1) < 0.15
)

# Dosing: 60 mg/m^2 over 1.5-hour IV infusion to central. amt in mg,
# rate in mg/h. Sampling: pre-dose, end-of-infusion, then the Methods
# schedule midpoints plus a 48-h tail for terminal-slope estimation.
infusion_h    <- 1.5
dose_per_bsa  <- 60                    # mg/m^2
sample_times  <- c(0, 1.5, 3, 6, 14, 24, 36, 48)

make_subject <- function(i) {
  row  <- cov_tab[i, ]
  amt  <- dose_per_bsa * row$BSA       # mg total
  rate <- amt / infusion_h             # mg/h
  doses <- tibble::tibble(
    id   = row$id,    time = 0,
    evid = 1L,        amt  = amt,
    rate = rate,      dv   = NA_real_
  )
  obs <- tibble::tibble(
    id   = row$id,    time = sample_times,
    evid = 0L,        amt  = NA_real_,
    rate = NA_real_,  dv   = NA_real_
  )
  bind_rows(doses, obs) |>
    mutate(
      WT  = row$WT,
      HT  = row$HT,
      BSA = row$BSA
    ) |>
    arrange(time, desc(evid))
}

events <- bind_rows(lapply(seq_len(n_subj), make_subject))

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

# Simulation

```{r simulate}
mod         <- readModelDb("Hempel_2003_daunorubicin_liposomal")
mod_typical <- rxode2::zeroRe(mod)

sim_typical <- rxode2::rxSolve(
  object = mod_typical, events = events,
  keep   = c("WT", "HT", "BSA")
) |>
  as.data.frame()

sim_stoch <- rxode2::rxSolve(
  object = mod, events = events,
  keep   = c("WT", "HT", "BSA")
) |>
  as.data.frame()
```

# Replicate published figures

## Figure 5 -- representative subject concentration-time profile

Hempel 2003 Figure 5 shows the measured plasma concentrations and
individual model predictions (model 15) for one representative patient
(ID 23) sampled across four administrations of Daunoxome. Plasma
concentrations during the 1.5-h infusion rise to a peak of order ~10000
ng/mL at end of infusion, then decline mono-exponentially with an
apparent terminal half-life of order 5-10 h. The figure below shows the
typical-value (`zeroRe`) prediction for the median 48.8 kg, 1.48 m^2
virtual subject, alongside the cohort-wide stochastic spread.

```{r figure-5}
typical_subject_id <- cov_tab |>
  arrange(abs(WT - 48.8)) |>
  slice(1) |>
  pull(id)

sim_typical_one <- sim_typical |>
  filter(id == typical_subject_id, time > 0)

ggplot() +
  geom_line(
    data = sim_stoch |> filter(time > 0),
    aes(time, Cc, group = id),
    colour = "gray70", alpha = 0.4
  ) +
  geom_line(
    data = sim_typical_one,
    aes(time, Cc),
    colour = "#1f77b4", linewidth = 1.1
  ) +
  geom_vline(xintercept = infusion_h, linetype = "dashed",
             colour = "gray40") +
  scale_y_log10() +
  labs(
    x = "Time after start of infusion (h)",
    y = "Total daunorubicin (ng/mL, log scale)",
    title = "Figure 5 (analogue) -- single-dose Daunoxome 60 mg/m^2",
    subtitle = paste(
      "Replicates Hempel 2003 Figure 5 schematic for a representative",
      "subject. Blue: typical-value (zeroRe) prediction for the",
      "median-weight (~48.8 kg) virtual subject. Gray:",
      "between-subject stochastic spread. Dashed line marks end of",
      "the 1.5-h infusion."
    )
  ) +
  theme_minimal()
```

## Concentration-time profile by body-weight tertile

The model parameterises CL and V as a linear product of body weight (no
reference-weight normalisation), so smaller patients have lower CL and V
in absolute terms and the per-dose AUC and t1/2 are roughly invariant of
weight in mass-normalised units. The plot below stratifies the
virtual cohort by weight tertile to show this directly.

```{r profile-by-weight}
wt_breaks <- quantile(cov_tab$WT, probs = c(0, 1/3, 2/3, 1))
sim_stoch <- sim_stoch |>
  mutate(wt_tert = cut(WT, breaks = wt_breaks,
                       include.lowest = TRUE,
                       labels = c("Low (lightest 1/3)",
                                  "Mid",
                                  "High (heaviest 1/3)")))

sim_stoch |>
  filter(time > 0) |>
  group_by(time, wt_tert) |>
  summarise(
    q05 = quantile(Cc, 0.05, na.rm = TRUE),
    q50 = quantile(Cc, 0.50, na.rm = TRUE),
    q95 = quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  ggplot(aes(time, q50, fill = wt_tert)) +
  geom_ribbon(aes(ymin = q05, ymax = q95), alpha = 0.2) +
  geom_line(aes(colour = wt_tert), linewidth = 0.9) +
  geom_vline(xintercept = infusion_h, linetype = "dashed",
             colour = "gray40") +
  scale_y_log10() +
  labs(
    x = "Time after start of infusion (h)",
    y = "Total daunorubicin (ng/mL, log scale)",
    colour = "Weight tertile",
    fill   = "Weight tertile",
    title  = "Concentration-time profile by body-weight tertile",
    subtitle = paste(
      "Median (line) and 5th-95th percentile (ribbon) of cohort",
      "stochastic simulation at 60 mg/m^2 IV infusion over 1.5 h."
    )
  ) +
  theme_minimal()
```

# PKNCA validation

The Hempel 2003 paper does not tabulate NCA parameters per subject;
however, Table 3 compares this study's BSA-normalised V, CL, terminal
half-life, and AUC at 60 mg/m^2 against four prior Daunoxome
publications. PKNCA is used here to summarise Cmax, AUCinf, half-life,
and clearance from the stochastic simulation, and the medians are
compared against the Table 3 "this study" column.

```{r pknca}
sim_for_nca <- sim_stoch |>
  filter(!is.na(Cc), Cc > 0, time > 0) |>
  select(id, time, Cc) |>
  mutate(treatment = "60 mg/m^2") |>
  as.data.frame()

doses_for_nca <- events |>
  filter(evid == 1L) |>
  select(id, time, amt) |>
  mutate(treatment = "60 mg/m^2") |>
  as.data.frame()

conc_obj <- PKNCA::PKNCAconc(
  data    = sim_for_nca,
  formula = Cc ~ time | treatment + id,
  concu   = "ng/mL",
  timeu   = "hr"
)
dose_obj <- PKNCA::PKNCAdose(
  data    = doses_for_nca,
  formula = amt ~ time | treatment + id,
  doseu   = "mg"
)

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

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

knitr::kable(
  summary(nca_res),
  caption = paste(
    "Simulated NCA parameters (PKNCA) for a single 60 mg/m^2 Daunoxome",
    "IV infusion over 1.5 h in the 24-subject virtual cohort."
  )
)
```

## Comparison against Table 3 (this-study column)

Hempel 2003 Table 3 reports, for the 60 mg/m^2 single-dose case
(BSA-normalised):

* V       = 1.93 L/m^2
* CL      = 0.233 L/h/m^2
* t1/2    = 5.66 h
* AUC at 60 mg/m^2 = 231.3 mg.h/L (= 231300 ng.h/mL)

The model parameterises CL and V on body weight, not BSA. For the
median paediatric subject (48.8 kg, 1.48 m^2, dose 88.8 mg), the
typical-value (`zeroRe`) prediction is

* CL_typ  = 0.00641 L/h/kg * 48.8 kg = 0.313 L/h, equivalent to
  0.313 / 1.48 = 0.211 L/h/m^2 (about 9% below Table 3's 0.233 L/h/m^2)
* V_typ   = 0.0654 L/kg * 48.8 kg = 3.19 L, equivalent to
  3.19 / 1.48 = 2.16 L/m^2 (about 12% above Table 3's 1.93 L/m^2)
* t1/2    = ln(2) * V_typ / CL_typ = 7.07 h (about 25% above Table
  3's 5.66 h)
* AUC     = 88.8 mg / 0.313 L/h = 283.5 mg.h/L (about 23% above
  Table 3's 231.3 mg.h/L)

The differences between the model-15 (weight-based, this packaged
implementation) and Table 3 (BSA-normalised individual NCA values for
the same cohort) reflect the different parameterisations the source
paper used for the two analyses: the per-BSA values in Table 3 were
computed for comparison against the prior literature (which also used
BSA-normalised reporting), whereas Table 2 model 15 is the per-weight
final population model that achieved the lowest objective function
(OF = 417). Both numbers are correct in their respective frames; the
packaged model returns the weight-based model-15 predictions.

```{r table-3-comparison}
typ_wt   <- 48.8
typ_bsa  <- 1.48
typ_dose <- 60 * typ_bsa
typ_cl_lh   <- 0.00641 * typ_wt
typ_v_l     <- 0.0654  * typ_wt
typ_t12_h   <- log(2)  * typ_v_l / typ_cl_lh
typ_auc_mghl <- typ_dose / typ_cl_lh
typ_cl_per_bsa <- typ_cl_lh / typ_bsa
typ_v_per_bsa  <- typ_v_l   / typ_bsa

published <- tibble::tibble(
  parameter         = c("V (L/m^2)", "CL (L/h/m^2)",
                        "t1/2 (h)", "AUC at 60 mg/m^2 (mg.h/L)"),
  table_3_value     = c(1.93, 0.233, 5.66, 231.3),
  model_15_typical  = c(round(typ_v_per_bsa, 3),
                        round(typ_cl_per_bsa, 3),
                        round(typ_t12_h, 2),
                        round(typ_auc_mghl, 1))
) |>
  mutate(pct_diff = round(100 * (model_15_typical / table_3_value - 1), 1))

knitr::kable(
  published,
  caption = paste(
    "Hempel 2003 Table 3 (this-study column, BSA-normalised) vs",
    "packaged model 15 typical-value prediction at the median",
    "paediatric subject (48.8 kg, 1.48 m^2). Differences",
    "reflect the per-weight vs per-BSA parameterisation used for",
    "the two analyses; the packaged model implements model 15",
    "(the paper's chosen final population model)."
  )
)
```

# Assumptions and deviations

- **Inter-occasion variability (IOV) of 16.7% CV on CL is NOT encoded
  structurally.** Hempel 2003 Table 2 model 15 retains an IOV
  component on CL alongside the IIV on CL (51% CV) and IIV on V (27%
  CV). The IOV captures within-subject variability in apparent
  clearance from one dosing occasion to the next (within a multi-cycle
  treatment course). nlmixr2lib has no canonical occasion-column
  convention, and the prevailing precedent for IOV without an
  operational OCC column is to omit it from the packaged model
  (Andrews 2017 tacrolimus, Brooks 2021 tacrolimus). Downstream users
  who want to reproduce the IOV can add an OCC indicator and a
  per-occasion eta in rxode2; the IOV variance on the log scale is
  log(1 + 0.167^2) = 0.0276. Population-level summaries (typical Cmax,
  typical AUC) are unaffected; only within-subject variability across
  cycles is missing.

- **Linear weight scaling with no reference-weight normalisation
  (per the source paper).** Model 15 uses CL = theta_CL * WT and V =
  theta_V * WT directly, with theta_CL = 6.41 mL/h/kg and theta_V =
  65.4 mL/kg reported per kg body weight. This is equivalent to a
  fixed allometric exponent of 1 with no reference weight; downstream
  comparison against models that use the more common allometric form
  (WT/70)^0.75 must account for the different functional shape.

- **Table 3 versus Table 2 model 15.** The paper reports two
  parameterisations of the cohort: Table 2 model 15 (weight-based,
  the chosen final model, OF = 417) and Table 3 (per-BSA values
  calculated specifically for comparison against the prior literature
  that also reported per-BSA values). The packaged model implements
  Table 2 model 15; the Table 3 comparison values are approximately
  20-25% different in t1/2 and AUC for the median paediatric subject
  because the per-BSA Table 3 values were computed for a different
  reporting frame, not because either set of numbers is wrong. See
  the comparison table above for the side-by-side breakdown.

- **Sex distribution not encoded.** Hempel 2003 reports that boys
  were heavier on average (median 58.0 vs 37.5 kg) but found no
  significant sex difference in PK after accounting for weight; the
  paper does not tabulate the exact M/F count. The
  `population$sex_female_pct` field is left as `NA` and sex is not a
  covariate of the model.

- **Body weight time-fixed in the simulation.** The paper does not
  explicitly state whether body weight was updated between dosing
  cycles. In the virtual cohort above weight is held constant across
  the 48-h observation window (a reasonable approximation given the
  short window). For multi-cycle simulations over weeks-to-months,
  users should consider supplying time-varying weight.

- **BSA computed via the Mosteller formula in the virtual cohort.**
  The source paper does not state which BSA formula it used; the
  Mosteller formula (BSA = sqrt(WT * HT / 3600), HT in cm, WT in kg)
  is applied here to derive BSA from the simulated WT and HT
  distributions for the 60 mg/m^2 dose calculation. The model itself
  does NOT use BSA as a covariate -- BSA appears only in the virtual
  cohort's dose-amount computation.

- **Concentration units (ng/mL) require an explicit *1000 scaling
  in `model()`.** With dose in mg and `vc` in L, the ratio
  `central / vc` carries units of mg/L = 1000 ng/mL. The `Cc` output
  multiplies by 1000 to express concentrations in the paper's
  reported ug/L units (numerically identical to ng/mL).
  `checkModelConventions()` may issue an info-level message flagging
  this dosing-vs-concentration magnitude mismatch; the scaling is
  intentional and documented in the model file.

- **Single-dose simulation only.** This vignette simulates one
  cycle's day-1 dose. The paper's dosing schedule (Daunoxome on days
  1 and 5 per cycle, four patients escalated to days 1, 3, 5) is
  mechanically straightforward to extend by adding further dose rows;
  the single-dose case is shown here as the validation pipeline.
