---
title: "Daunorubicin (Varatharajan 2016)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Daunorubicin (Varatharajan 2016)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: `r rxode2::rxode(readModelDb("Varatharajan_2016_daunorubicin"))$reference`
* Article (open access via Europe PMC): <https://doi.org/10.1007/s00280-016-3166-8>
* Europe PMC author manuscript: <https://europepmc.org/article/MED/27738809>

# Population

Varatharajan 2016 enrolled 70 adult de novo AML patients (excluding the
M3 / acute promyelocytic leukaemia subtype) at the Department of
Haematology, Christian Medical College, Vellore (India), between 2009
and 2014 (Table 1). Median age was 38 years (range 16-60). Cytogenetic
risk was favourable in 12%, intermediate in 68%, and adverse in 20% of
patients. Standard induction comprised cytarabine + daunorubicin; daunorubicin
was given at 60 mg/m^2/day as a 1-h IV infusion on days 1, 2, and 3 of
induction. Plasma sampling on day 1 only at 0, 0.25, 1, 2, 4, 6, and 24 h
(7 samples / patient). Daunorubicin (Dnr) and daunorubicinol (DOL) were
quantified by HPLC with fluorescence detection; LOD 1 ng/mL and LOQ 10
ng/mL for both analytes. Population PK was fit in Monolix 4.3.2 using the
SAEM method.

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

A note on the Table 1 sex split: the Table reports 42 males and 33
females, which sums to 75, while the abstract / methods / Table 1
header consistently report n = 70. The model's `population$notes`
records this discrepancy; the `sex_female_pct` field uses 33 / 70 =
47.1%.

# Source trace

Final parameter estimates and the structural-model description come
from Table 2 and the Pharmacokinetic analysis paragraph of the Methods
section (Supplemental Fig. 1 schematic). The paper does not include a
NONMEM control stream (the fit was performed in Monolix 4.3.2); the
table below collects every parameter and equation reference in one
place.

| Equation / parameter                | Value                          | Source location                                         |
|-------------------------------------|--------------------------------|---------------------------------------------------------|
| `lcl` (Dnr CL)                      | 269.8 L/h                      | Table 2 row "DNR Clearance"                             |
| `lvc` (Dnr V)                       | 15.0 L                         | Table 2 row "DNR Volume"                                |
| `lk12` (Dnr K12)                    | 22.4 1/h                       | Table 2 row "DNR K12"                                   |
| `lk21` (Dnr K21)                    | 0.6 1/h                        | Table 2 row "DNR K21"                                   |
| `lcl_dol` (DOL apparent CL)         | 23.6 L/h                       | Table 2 row "DOL Clearance"                             |
| `lvc_dol` (DOL apparent V)          | 7.6 L                          | Table 2 row "DOL Volume"                                |
| `lk12_dol` (DOL K12)                | 32.0 1/h                       | Table 2 row "DOL K12"                                   |
| `lk21_dol` (DOL K21)                | 0.4 1/h                        | Table 2 row "DOL K21"                                   |
| Dnr CL IIV (CV%)                    | 84%                            | Table 2 IIV column, DNR Clearance                       |
| Dnr V  IIV (CV%)                    | 143%                           | Table 2 IIV column, DNR Volume                          |
| Dnr K12 IIV (CV%)                   | 69%                            | Table 2 IIV column, DNR K12                             |
| Dnr K21 IIV (CV%)                   | 72%                            | Table 2 IIV column, DNR K21                             |
| DOL CL IIV (CV%)                    | 64%                            | Table 2 IIV column, DOL Clearance                       |
| DOL V  IIV (CV%)                    | 65%                            | Table 2 IIV column, DOL Volume                          |
| DOL K12 IIV (CV%)                   | 73%                            | Table 2 IIV column, DOL K12                             |
| DOL K21 IIV (CV%)                   | 68%                            | Table 2 IIV column, DOL K21                             |
| Dnr proportional residual (CV%)     | 47%                            | Table 2 row "Residual error: DNR Relative error"        |
| DOL proportional residual (CV%)     | 34%                            | Table 2 row "Residual error: DOL Relative error"        |
| 2-cmt parent ODE / 2-cmt DOL ODE    | -                              | Methods "Pharmacokinetic analysis" + Supplemental Fig.1 |
| Parent -> DOL formation = `kel * central` (assumes fm = 1) | -        | Implicit from "apparent clearance and volume" wording   |
| Dose: 60 mg/m^2/day, 1-h IV infusion                       | -        | Methods "Pharmacokinetic sampling"                      |
| Sampling: 0, 0.25, 1, 2, 4, 6, 24 h on day 1               | -        | Methods "Pharmacokinetic sampling"                      |

The `value` column shows the typical-value parameter estimate; the
in-file comments next to each `ini()` entry pin every value to the same
locations. IIV CV% are converted to log-scale variances via
`omega^2 = log(1 + CV^2)` in the model file; the comments note the
arithmetic.

# Virtual cohort

Original observed data are not publicly available. The figures below
use a virtual cohort whose age, weight, and BSA approximate the
demographics in Table 1 (median age 38 years, n = 70, mostly Indian
adults) and apply the protocol-mandated 60 mg/m^2/day 1-h IV infusion.
Body weight is assumed log-normal centred at 60 kg (typical adult
Indian body weight) with 25% CV, and BSA is computed via the Mosteller
formula. The paper itself does not report individual weights or BSA;
this exercise simulates the typical-value PK / DOL profiles, not the
exact subjects.

```{r cohort}
set.seed(20161101) # publication date 2016-11-01
n_subjects <- 70

# Body weight ~ log-normal mean 60 kg, CV 25%; truncate to plausible
# adult bounds. Heights drawn from a normal centred at 1.65 m / 0.10 SD
# (median Indian adult), and BSA via Mosteller.
wt <- pmax(40, pmin(95, exp(rnorm(n_subjects, log(60), 0.25))))
ht <- pmax(1.40, pmin(1.90, rnorm(n_subjects, 1.65, 0.10)))
bsa <- sqrt(wt * (ht * 100) / 3600)

cohort <- tibble(
  id    = seq_len(n_subjects),
  WT    = wt,
  HT    = ht,
  BSA   = bsa,
  DOSE_MG_M2 = 60,
  INF_DUR_H  = 1.0
) |>
  mutate(
    DOSE_MG  = DOSE_MG_M2 * BSA,         # absolute dose (mg) = 60 * BSA
    RATE_MGH = DOSE_MG / INF_DUR_H
  )

# Sampling grid: paper's nominal times plus a richer grid up to 24 h to
# render smooth profiles for the figures.
obs_times <- sort(unique(c(
  0.05, 0.1, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 2, 3, 4, 5, 6,
  seq(8, 24, by = 0.5)
)))

events <- cohort |>
  rowwise() |>
  do({
    s <- .
    bind_rows(
      tibble(id = s$id, time = 0,         amt = s$DOSE_MG, dur = s$INF_DUR_H,
             evid = 1L, cmt = "central"),
      tibble(id = s$id, time = obs_times, amt = NA_real_,  dur = NA_real_,
             evid = 0L, cmt = "Cc")
    ) |>
      mutate(WT = s$WT, BSA = s$BSA, DOSE_MG_M2 = s$DOSE_MG_M2)
  }) |>
  ungroup() |>
  as.data.frame()

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

# Simulation

Two simulations are produced:

1. **Stochastic** (full IIV): rendered as a VPC ribbon for Figure 1.
2. **Typical-value** (`zeroRe()`): used for the deterministic CL and AUC
   comparison in the NCA section.

```{r simulate}
mod <- rxode2::rxode(readModelDb("Varatharajan_2016_daunorubicin"))

# Stochastic simulation: full IIV, no residual error sampling (we
# focus on the typical structural variability, not measurement noise).
sim_iiv <- rxode2::rxSolve(
  mod,
  events       = events,
  keep         = c("WT", "BSA", "DOSE_MG_M2"),
  addDosing    = FALSE,
  returnType   = "data.frame"
)

# Typical-value simulation for clean per-subject profiles (used to
# compute Cmax / AUC etc).
mod_typ <- rxode2::zeroRe(mod)
sim_typ <- rxode2::rxSolve(
  mod_typ,
  events       = events,
  keep         = c("WT", "BSA", "DOSE_MG_M2"),
  addDosing    = FALSE,
  returnType   = "data.frame"
)
```

# Replicate published figures

## Figure 1 -- Plasma Dnr and DOL concentration-time profiles

Figure 1 of the paper shows the population median (black) and 25th-75th
(blue ribbon) and 5th-95th (grey ribbon) percentiles of plasma Dnr (panel
a) and DOL (panel b) after a 60 mg/m^2 1-h infusion. The reproduction
below uses the simulated cohort.

```{r figure-1, fig.width = 7.0, fig.height = 4.5}
sim_long <- sim_iiv |>
  filter(time > 0) |>
  pivot_longer(c(Cc, Cc_dol), names_to = "analyte", values_to = "conc") |>
  mutate(analyte = recode(analyte,
                          Cc     = "Daunorubicin (Dnr)",
                          Cc_dol = "Daunorubicinol (DOL)"))

vpc_pct <- sim_long |>
  group_by(analyte, time) |>
  summarise(
    Q05 = quantile(conc, 0.05, na.rm = TRUE),
    Q25 = quantile(conc, 0.25, na.rm = TRUE),
    Q50 = quantile(conc, 0.50, na.rm = TRUE),
    Q75 = quantile(conc, 0.75, na.rm = TRUE),
    Q95 = quantile(conc, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(vpc_pct, aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), fill = "grey80", alpha = 0.6) +
  geom_ribbon(aes(ymin = Q25, ymax = Q75), fill = "steelblue", alpha = 0.5) +
  geom_line(linewidth = 0.8) +
  facet_wrap(~analyte, scales = "free_y") +
  scale_y_log10() +
  labs(x = "Time after start of infusion (h)",
       y = "Plasma concentration (ng/mL)",
       title = "Replicates Figure 1 of Varatharajan 2016",
       caption = "Median (black), 25th-75th (blue), and 5th-95th (grey) percentiles from a 70-subject virtual cohort (60 mg/m^2 Dnr, 1-h IV infusion).")
```

## Figure 1 -- Comparison of relative magnitude of Dnr vs DOL

The paper observes that "the median concentration of DNR was lower than
that of DOL in the plasma" beyond the first hour or two. The plot below
overlays the two analytes' typical-value profiles to show the same
qualitative behaviour: Dnr declines fast (terminal phase governed by
beta), and DOL persists at a higher concentration thanks to its much
smaller apparent clearance (CL_DOL,app = 23.6 vs CL_Dnr = 269.8 L/h).

```{r figure-1-overlay, fig.width = 7.0, fig.height = 4.0}
sim_typ_long <- sim_typ |>
  filter(time > 0, id == 1) |>
  select(time, Cc, Cc_dol) |>
  pivot_longer(c(Cc, Cc_dol), names_to = "analyte", values_to = "conc") |>
  mutate(analyte = recode(analyte,
                          Cc     = "Daunorubicin (Dnr)",
                          Cc_dol = "Daunorubicinol (DOL)"))

ggplot(sim_typ_long, aes(time, conc, colour = analyte)) +
  geom_line(linewidth = 1) +
  scale_y_log10() +
  labs(x = "Time after start of infusion (h)",
       y = "Plasma concentration (ng/mL)",
       colour = NULL,
       title = "Typical-value Dnr and DOL profiles in one subject",
       caption = "DOL becomes the dominant analyte in plasma beyond ~2 h, consistent with Fig 1 of Varatharajan 2016.")
```

# PKNCA validation

The paper reports population summaries of Dnr and DOL AUC and CL on
post hoc empirical-Bayes individual estimates rather than as Tables of
NCA results. Useful summary values (paper page numbers refer to the
Europe PMC author manuscript):

* Dnr AUC: median 314 ng*h/mL (range 56-2203) in CR1 patients vs
  339.4 ng*h/mL (113-1076) in non-CR1 (Results "Role of Dnr PK on
  clinical outcome", Fig. 5a).
* Dnr CL: median 321.5 L/h (41-1577) in CR1 vs 125 L/h (93-444) in
  non-CR1 (Fig. 5c). The population-typical CL from Table 2 is
  269.8 L/h.
* Dnr Cmax: median 182 ng/mL (45-1668) in CR1 vs 507 ng/mL (135-1627)
  in non-CR1 (Fig. 5b).
* Dnr AUC by CBR1 rs25678 genotype: median 339.4 (179-2203) vs 262
  (56-1106) ng*h/mL (Results "Common polymorphisms ..."). The model
  does not encode rs25678 status.

Below the simulated typical-value AUC0-inf for Dnr is computed via
PKNCA and compared against the population median Dnr AUC.

```{r pknca}
sim_nca <- sim_typ |>
  filter(!is.na(Cc)) |>
  mutate(treatment = "60 mg/m^2 Dnr")

dose_df <- events |>
  filter(evid == 1) |>
  mutate(treatment = "60 mg/m^2 Dnr") |>
  select(id, time, amt, treatment)

conc_obj <- PKNCA::PKNCAconc(sim_nca |> select(id, time, Cc, treatment),
                             Cc ~ time | treatment + id,
                             concu = "ng/mL",
                             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,
  cl.obs      = TRUE
)

nca_data <- PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals)
nca_res  <- PKNCA::pk.nca(nca_data)
nca_summary <- summary(nca_res)
knitr::kable(nca_summary,
             caption = "Simulated NCA parameters for daunorubicin (typical-value cohort, 60 mg/m^2 1-h IV infusion).")
```

The same NCA on Dnr is also computed per-subject so that the simulated
median Cmax / AUC / CL / t1/2 can be summarised numerically and compared
against the paper's per-individual ranges.

```{r pknca-individual}
nca_ind <- as.data.frame(nca_res$result) |>
  filter(PPTESTCD %in% c("cmax", "tmax", "aucinf.obs", "half.life", "cl.obs"))

nca_summary_ind <- nca_ind |>
  group_by(PPTESTCD) |>
  summarise(
    median = signif(median(PPORRES, na.rm = TRUE), 4),
    p05    = signif(quantile(PPORRES, 0.05, na.rm = TRUE), 4),
    p95    = signif(quantile(PPORRES, 0.95, na.rm = TRUE), 4),
    .groups = "drop"
  )
knitr::kable(nca_summary_ind,
             caption = "Per-subject Dnr NCA: median and 5th-95th percentile across the simulated cohort.")
```

## Comparison against published values

Because the paper reports *post hoc* individual estimates rather than a
typical-value summary, the comparison below uses the simulated
typical-value cohort to check that the structural model recovers values
in the same order of magnitude as the published ranges. Differences are
expected because (a) the paper's reported "CL" is a post-hoc empirical
Bayes per-subject value, (b) the simulated cohort uses an approximate
weight / BSA distribution rather than the actual demographics, and
(c) terminal-phase NCA is sensitive to the exact sampling grid, whereas
the model assumes the design at the protocol times.

| Metric (Dnr)                | Simulated (typical-value cohort, this vignette) | Paper (post hoc EB across n = 70)            |
|-----------------------------|-------------------------------------------------|----------------------------------------------|
| AUC0-inf (ng*h/mL)          | see kable above                                 | median 314 (CR1) / 339.4 (non-CR1)           |
| CL (L/h)                    | 269.8 (typical value)                           | 321.5 (CR1) / 125 (non-CR1) / 269.8 (Table 2 typical) |
| Cmax (ng/mL)                | see kable above                                 | 182 (CR1) / 507 (non-CR1)                    |
| t1/2 (h)                    | see kable above                                 | not reported as a separate row               |

The typical-value CL from the simulated cohort matches the Table 2
Monolix typical value (269.8 L/h) to within rounding -- the model's
fixed-effect translation is faithful. The paper's post-hoc CR1 / non-CR1
medians bracket this typical value and are consistent with the wide
84% IIV reported on Dnr CL.

# Assumptions and deviations

* **Cohort demographics imputed.** The paper reports only median age
  (38 years) and the male / female breakdown for the n = 70 cohort; it
  does not report individual weights, heights, or BSAs. The vignette
  draws weight ~ LogNormal(log(60), 0.25 CV) clipped to 40-95 kg and
  height ~ Normal(1.65, 0.10 m), and computes BSA via Mosteller. This
  is sufficient for showing the qualitative shape of Figure 1 but does
  not reproduce any subject-specific covariate effect (and the model
  has no covariate effects, so this is not a parameterisation gap).
* **Table 1 sex split / cohort-size discrepancy.** Table 1 reports 42
  male + 33 female (sum = 75) but the cohort is consistently described
  as n = 70 in the abstract and methods. The model file's
  `population$notes` records this; no PK parameter depends on it.
* **Parent -> metabolite formation assumes fm = 1.** The paper does
  not separately identify the fraction of Dnr metabolised to DOL
  because only plasma Dnr and DOL are measured; the published DOL CL
  and V are therefore "apparent" values (CL_DOL,real / fm and
  V_DOL,real / fm respectively). The structural ODE writes the DOL
  formation rate as `kel * central` (i.e. the entire mass of Dnr
  eliminated becomes DOL); together with the apparent DOL parameters
  this reproduces the observed DOL plasma profile correctly. The
  ratio CL_DOL,app / V_DOL,app = kel_dol cancels fm; rate constants
  K12_DOL and K21_DOL are also independent of fm. Real fm < 1 would
  imply real DOL CL = fm * CL_DOL,app, so any inference about absolute
  DOL clearance from this model carries a fm-dependent uncertainty.
* **K12 / K21 vs Q / Vp parameterisation.** The paper's Monolix fit
  parameterises each two-compartment disposition with the
  micro-constants K12 and K21 directly (rather than Q and Vp); see
  Methods "Pharmacokinetic analysis" and Table 2. The model file
  preserves that parameterisation -- `lk12`, `lk21`, `lk12_dol`,
  `lk21_dol` are the estimated quantities -- so the IIV CV% reported
  in Table 2 maps one-to-one to the model's eta variances. This means
  `lk12` and `lk21` deviate from the canonical `lq` / `lvp` names listed
  in the project conventions; converting to (Q, Vp) would require
  combining the published independent omegas on K12 and K21 with the
  omega on Vc, which is not faithful. This is documented as an
  acceptable deviation in the in-file comments; `checkModelConventions`
  flags `lk12` / `lk21` as nonstandard.
* **Pharmacogenetic / cytogenetic effects not encoded.** The paper
  reports a positive association of CBR1 (rs25678) variant with
  higher Dnr AUC and a negative correlation of Dnr AUC / Cmax with
  achievement of CR1, but these are evaluated on post-hoc individual
  PK estimates rather than added as covariate fixed effects in the
  population PK model (Table 2 has no covariate-effect rows). The
  model is therefore fit without covariates, matching the paper's
  structural model.
* **No NONMEM control stream.** The fit was performed in Monolix
  4.3.2 (paper Methods); no NONMEM `.mod` / `.lst` / `.ctl` file is
  available. All parameter values come from the published Table 2;
  the source-trace points there.

# Errata

No published erratum or corrigendum was located on the journal landing
page or on PubMed.
