---
title: "Eribulin mesilate (van Hasselt 2013)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Eribulin mesilate (van Hasselt 2013)}
  %\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

* Citation: `r readModelDb("vanHasselt_2013_eribulin")()$meta$reference`
* Description: `r readModelDb("vanHasselt_2013_eribulin")()$meta$description`
* Article: <https://doi.org/10.1111/bcp.12143>

The model file at `inst/modeldb/specificDrugs/vanHasselt_2013_eribulin.R` couples
the Majid 2014 three-compartment eribulin popPK (fixed) with the van Hasselt 2013
Friberg-style semi-physiological PD model for chemotherapy-induced neutropenia.
The PK layer's parameter values are FIXED from Majid 2014 because van Hasselt
2013 explicitly cites the same Eisai popPK analysis as "unpublished data on
file, Dr Z. Hussein" and does not report PK parameter values in the paper (see
the Errata section for the substitution rationale). The PD layer's parameters
come from van Hasselt 2013 Table 3 "Full covariate model" column.

# Population

The van Hasselt 2013 pooled data set combines 12 phase I / II / III studies of
eribulin mesilate for a total of 1579 patients with late-stage metastatic breast
cancer, contributing 23 427 absolute neutrophil count (ANC) measurements. PK
data were available for 428 patients (27%). The cohort is predominantly female
(85.6%; 229 male / 1359 female per Table 2) and ethnically diverse: 79.4%
Caucasian, 6.1% Japanese, 5.3% Black/African American, 1.5% Asian/Pacific
Islander, 8.4% Other/unknown. Median age 58.0 years (IQR 49-66), median
bodyweight 67.7 kg (IQR 59-77.6), median BSA 1.73 m^2 (IQR 1.61-1.86).

Baseline laboratory covariates (Table 2 medians): albumin 3.90 g/dL,
alkaline phosphatase 118 U/L, alanine transaminase 25 U/L, aspartate
transaminase 30 U/L, bilirubin 0.50 mg/dL, lactate dehydrogenase 328 U/L,
platelets 260 x 10^9/L, protein 7.1 g/dL. G-CSF was administered as concomitant
medication in 382 / 1302 (yes/no) patients (22.7% of the cohort). Prior therapy
exposure was high across the cohort: 96% previous chemotherapy, 28% previous
Pt-containing chemotherapy, 76% previous radiotherapy, 58% previous hormonal
therapy, 10% received blood transfusions.

The same information is available programmatically via
`readModelDb("vanHasselt_2013_eribulin")()$meta$population`.

# Source trace

The per-parameter origin is recorded as an in-file comment next to each `ini()`
entry. The table below collects the key entries in one place; parameter names
match those in `inst/modeldb/specificDrugs/vanHasselt_2013_eribulin.R`.

| Equation / parameter                    | Value          | Source location                                    |
|-----------------------------------------|----------------|----------------------------------------------------|
| PK CL (typical, L/h)                    | 3.11           | Majid 2014 (via Kawamura 2018 section 2.3)         |
| PK V1 (L)                               | 4.06           | Majid 2014                                         |
| PK Q2, V2, Q3, V3                       | 2.64, 2.42, 6.60, 121 | Majid 2014                                  |
| PK covariate exponents (WT/ALB/ALP/TBILI on CL) | 0.75, 0.946, -0.209, -0.180 | Majid 2014                    |
| PK V1/V2/V3 linear WT scaling          | exponent 1     | Majid 2014                                         |
| PD ANC0 (baseline, cells/uL)            | 4030           | van Hasselt 2013 Table 3 final (4.03e9 cells/L)    |
| PD MTT (h)                              | 109            | van Hasselt 2013 Table 3 final                     |
| PD gamma (unitless)                     | 0.216          | van Hasselt 2013 Table 3 final                     |
| PD SLOPE (L/ug)                         | 0.0451         | van Hasselt 2013 Table 3 final                     |
| ALB power on MTT / SLOPE (ref 4 g/dL)   | 0.374 / 0.763  | van Hasselt 2013 Table 3 final                     |
| TBILI power on MTT (ref 2 mg/dL)        | -0.046         | van Hasselt 2013 Table 3 final                     |
| ALP power on MTT (ref 100 U/L)          | -0.0337        | van Hasselt 2013 Table 3 final                     |
| LDH power on MTT (ref 238 U/L)          | -0.0561        | van Hasselt 2013 Table 3 final                     |
| AST power on SLOPE (ref median 30 U/L)  | 0.119          | van Hasselt 2013 Table 3 final                     |
| CONMED_GCSF proportional on MTT         | 0.883          | van Hasselt 2013 Table 3 final                     |
| CONMED_GCSF proportional on SLOPE       | 1.3            | van Hasselt 2013 Table 3 final                     |
| IIV ANC0 / MTT / gamma / SLOPE (CV%)    | 37.3 / 13.9 / 12.2 / 54.0 | van Hasselt 2013 Table 3 final          |
| Proportional residual RUV (CV%)         | 49.6           | van Hasselt 2013 Table 3 final                     |
| PD equations (Prol + 3 transit + Circ)  | Eqs. 1-6       | van Hasselt 2013 Methods section                   |
| Drug-effect equation E = SLOPE * C      | Eq. 6          | van Hasselt 2013 Methods                           |
| Covariate model Eq. 10                  | product form   | van Hasselt 2013 Methods                           |

# Virtual cohort

The van Hasselt 2013 raw data are not publicly available. The virtual cohort
below approximates the Methods section 'Evaluation of the impact of covariates
on risk for neutropenia' simulation strategy: 200 patients (a smaller cohort
than the paper's 2000 to keep the vignette render under the 5-minute budget),
BSA sampled from N(mean = 1.57, sd = 0.22) truncated to [1, 3] m^2 per the
paper, weight fixed at the cohort median 67.7 kg, laboratory covariates set to
the Table 2 medians (ALB 3.90 g/dL = 39 g/L; ALP 118 U/L; TBILI 0.50 mg/dL =
8.55 umol/L; AST 30 U/L; LDH 328 U/L), and G-CSF-receival at 22.7% prevalence.
Approved eribulin dosing schedule: 1.4 mg/m^2 mesilate (= 1.23 mg/m^2 free base
via the Majid 2014 / Kawamura 2018 conversion factor 1.23/1.4) on day 1 and day
8 of a 21-day cycle.

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

# Truncated normal BSA per the paper simulation Methods
bsa_draw <- function(n, mean = 1.57, sd = 0.22, lo = 1, hi = 3) {
  x <- rnorm(2 * n, mean = mean, sd = sd)
  x <- x[x >= lo & x <= hi]
  x[seq_len(n)]
}

cohort <- tibble(
  id  = seq_len(n_sub),
  BSA = bsa_draw(n_sub),
  WT  = 67.7,          # cohort median (Table 2)
  ALB = 39,            # 3.90 g/dL = 39 g/L (canonical SI)
  ALP = 118,           # U/L
  TBILI = 8.55,        # 0.50 mg/dL = 8.55 umol/L (canonical SI)
  AST = 30,            # U/L
  LDH = 328,           # U/L
  CONMED_GCSF = as.integer(runif(n_sub) < 0.227)  # 22.7% G-CSF prevalence
) |>
  mutate(dose_mg = 1.23 * BSA)  # 1.23 mg/m^2 free base * BSA

# Dosing schedule: day 1 (t=0 h) and day 8 (t=168 h) IV over 2 min = 1/30 h
dose_rows <- cohort |>
  transmute(id, WT, ALB, ALP, TBILI, AST, LDH, CONMED_GCSF, amt = dose_mg,
            cmt = "central", evid = 1L, time = 0, dur = 1/30) |>
  bind_rows(
    cohort |> transmute(id, WT, ALB, ALP, TBILI, AST, LDH, CONMED_GCSF,
                        amt = dose_mg, cmt = "central", evid = 1L,
                        time = 168, dur = 1/30)
  )

# Observation grid: every 3 h for 21 days
obs_grid <- expand_grid(
  id = cohort$id,
  time = seq(0, 24 * 21, by = 3)
) |>
  left_join(cohort |> select(-BSA, -dose_mg), by = "id") |>
  mutate(amt = 0, cmt = "central", evid = 0L, dur = 0)

events <- bind_rows(dose_rows, obs_grid) |>
  arrange(id, time, desc(evid))
```

# Simulation

```{r simulate}
mod <- readModelDb("vanHasselt_2013_eribulin")
sim <- rxode2::rxSolve(mod, events = events,
                       keep = c("CONMED_GCSF", "BSA"),
                       addDosing = FALSE) |>
  as.data.frame()
```

For deterministic replication (typical-value trajectories), zero out random
effects:

```{r simulate-typical}
mod_typical <- readModelDb("vanHasselt_2013_eribulin") |> rxode2::zeroRe()
sim_typical <- rxode2::rxSolve(mod_typical, events = events, addDosing = FALSE) |>
  as.data.frame()
```

# Replicate published figures

## Figure 4 - ANC time profile for the typical patient

Reproduces the qualitative shape of van Hasselt 2013 Figure 4 for the typical
patient (all covariates at reference values, no G-CSF): ANC drops from the
4030 cells/uL baseline into a nadir 8-14 days after each dose, then recovers
via the transit-chain feedback loop.

```{r figure-4-typical}
typ_med <- sim_typical |>
  filter(id == 1) |>
  select(time, ANC, Cc)

ggplot(typ_med, aes(time / 24, ANC)) +
  geom_line(colour = "firebrick", linewidth = 0.8) +
  geom_hline(yintercept = 1000, linetype = "dashed", colour = "grey40") +
  geom_hline(yintercept = 500,  linetype = "dashed", colour = "grey40") +
  annotate("text", x = 0.5, y = 1050, hjust = 0, label = "Grade 3 (ANC < 1000)",
           colour = "grey30", size = 3) +
  annotate("text", x = 0.5, y = 550,  hjust = 0, label = "Grade 4 (ANC < 500)",
           colour = "grey30", size = 3) +
  geom_vline(xintercept = c(0, 7), linetype = "dotted", colour = "steelblue") +
  scale_x_continuous(breaks = seq(0, 21, by = 3)) +
  scale_y_continuous(limits = c(0, 5000)) +
  labs(x = "Time since first dose (days)",
       y = "ANC (cells/uL)",
       title = "Typical-patient ANC trajectory for approved eribulin dosing (day 1, day 8 of a 21-day cycle)",
       caption = "Replicates the qualitative shape of van Hasselt 2013 Figure 4 (continuous red line = typical parameter values).")
```

## Figure 2 - Visual predictive check across 200 virtual subjects

Reproduces the general shape of van Hasselt 2013 Figure 2 (VPC of ANC vs
time by study). Because the raw study-stratified data are not on disk, the
figure below is a single-panel VPC across the full virtual cohort at reference
covariate values with 22.7% G-CSF prevalence.

```{r figure-2-vpc}
vpc_df <- sim |>
  group_by(time) |>
  summarise(
    Q05 = quantile(ANC, 0.05, na.rm = TRUE),
    Q50 = quantile(ANC, 0.50, na.rm = TRUE),
    Q95 = quantile(ANC, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(vpc_df, aes(time / 24, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), fill = "steelblue", alpha = 0.25) +
  geom_line(colour = "steelblue", linewidth = 0.8) +
  geom_hline(yintercept = 1000, linetype = "dashed", colour = "grey40") +
  geom_hline(yintercept = 500,  linetype = "dashed", colour = "grey40") +
  scale_x_continuous(breaks = seq(0, 21, by = 3)) +
  labs(x = "Time since first dose (days)",
       y = "ANC (cells/uL)",
       title = "VPC across 200 virtual subjects at reference covariate values",
       caption = "Model-predicted median (line) and 90% prediction interval (shaded). Replicates the qualitative shape of van Hasselt 2013 Figure 2.")
```

## Table 5 base-model grade 3 / 4 incidences

Table 5 of van Hasselt 2013 reports the incidences of grade 3 neutropenia
(ANC < 1000 cells/uL for > 7 days) and grade 4 neutropenia (ANC < 500 cells/uL
for > 7 days) at the approved 1.4 mg/m^2 dose using n = 2000 stochastic
simulations: 26.6% grade 3, 16.85% grade 4. The paper's Table 5 is derived from
the base (no-covariate) model; the model here uses the final covariate model
with covariates fixed at reference values, so the two are not identical
parameterisations. The comparison below is nonetheless a useful check that
the model's neutropenia depth and duration are in the expected clinical range.

```{r table-5-incidence}
grade_3_thr <- 1000
grade_4_thr <- 500
window_h   <- 24 * 7  # > 7 days = > 168 h

# For each subject and each grade, compute the longest run of ANC below the
# threshold and flag the subject as having that grade if the run exceeds 7 days.
run_below <- function(times, values, threshold) {
  below <- values < threshold
  # rle on the below-threshold logical, weighted by time deltas
  if (!any(below)) return(0)
  dt <- c(diff(times), 0)
  rl <- rle(below)
  ends <- cumsum(rl$lengths)
  starts <- ends - rl$lengths + 1
  dur <- mapply(function(s, e) sum(dt[s:e]), starts, ends)
  # only true runs count
  max(dur[rl$values], 0)
}

incidence <- sim |>
  arrange(id, time) |>
  group_by(id) |>
  summarise(
    dur_g3 = run_below(time, ANC, grade_3_thr),
    dur_g4 = run_below(time, ANC, grade_4_thr),
    .groups = "drop"
  ) |>
  summarise(
    grade_3_pct = mean(dur_g3 > window_h) * 100,
    grade_4_pct = mean(dur_g4 > window_h) * 100
  )

sim_vs_paper <- tibble::tribble(
  ~metric,             ~simulated,               ~paper_base_model, ~note,
  "Grade 3 incidence (% > 7 days ANC < 1000 cells/uL)", incidence$grade_3_pct, 26.6,
    "van Hasselt 2013 Table 5 base-model row",
  "Grade 4 incidence (% > 7 days ANC < 500 cells/uL)",  incidence$grade_4_pct, 16.85,
    "van Hasselt 2013 Table 5 base-model row"
)

sim_vs_paper |>
  dplyr::rename(
    "Metric" = metric,
    "Simulated (n = 200, this vignette)" = simulated,
    "van Hasselt 2013 Table 5" = paper_base_model,
    "Note" = note
  ) |>
  knitr::kable(digits = 1, align = c("l", "r", "r", "l"),
               caption = "Grade 3 and grade 4 neutropenia incidence at 1.4 mg/m^2 dose over 21-day cycle 1. The vignette simulates the FINAL COVARIATE model at reference covariate values (n = 200); the paper's Table 5 uses the BASE (no-covariate) model with n = 2000 and covariates absorbed into wider IIVs. Grade 3 matches within one sample-SE; grade 4 is markedly lower because the reference-covariate cohort under-represents the deep-nadir tail (see Assumptions).")
```

# PKNCA validation

The paper does not tabulate NCA parameters for eribulin plasma
concentration, so no side-by-side simulated-vs-published NCA table is
possible for this extraction. The PKNCA block below is still included as a
sanity check that the Majid 2014 PK layer produces realistic exposure
values (published eribulin Cmax at 1.4 mg/m^2 is ~150-200 ng/mL; AUC over
one dosing cycle is on the order of 500-1000 ng*h/mL depending on covariate
mix).

```{r pknca}
sim_nca <- sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc) |>
  dplyr::mutate(treatment = "1.4 mg/m^2")

# Time-zero row (Cc = 0 pre-dose)
sim_nca <- dplyr::bind_rows(
  sim_nca,
  sim_nca |> dplyr::distinct(id, treatment) |>
    dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, treatment, time, .keep_all = TRUE) |>
  dplyr::arrange(id, treatment, time)

dose_df <- dose_rows |>
  filter(time == 0) |>
  select(id, time, amt) |>
  mutate(treatment = "1.4 mg/m^2")

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

# AUC over one dosing interval (0-168 h = day 1 to day 8)
intervals <- data.frame(
  start = 0,
  end   = 168,
  cmax  = TRUE,
  tmax  = TRUE,
  auclast = TRUE,
  half.life = TRUE
)

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

nca_summary <- as.data.frame(nca_res$result) |>
  dplyr::filter(PPTESTCD %in% c("cmax", "tmax", "auclast", "half.life")) |>
  dplyr::group_by(PPTESTCD) |>
  dplyr::summarise(
    median = median(PPORRES, na.rm = TRUE),
    p05    = quantile(PPORRES, 0.05, na.rm = TRUE),
    p95    = quantile(PPORRES, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  dplyr::mutate(units = dplyr::case_when(
    PPTESTCD == "cmax"     ~ "mg/L (1 mg/L = 1000 ng/mL)",
    PPTESTCD == "tmax"     ~ "hour",
    PPTESTCD == "auclast"  ~ "mg*h/L",
    PPTESTCD == "half.life" ~ "hour",
    TRUE ~ ""
  ))

nca_summary |>
  dplyr::rename(
    "NCA parameter"                   = PPTESTCD,
    "Median (n = 200)"                = median,
    "5th percentile"                  = p05,
    "95th percentile"                 = p95,
    "Units"                           = units
  ) |>
  knitr::kable(digits = 4, align = c("l", "r", "r", "r", "l"),
               caption = "Simulated eribulin plasma NCA over dose interval 0-168 h at 1.4 mg/m^2 mesilate (= 1.23 mg/m^2 free base). Reported for sanity: expected clinical Cmax ~ 0.15-0.2 mg/L (150-200 ng/mL). The paper itself does not report per-subject NCA.")
```

# Assumptions and deviations

- **PK layer parameter values FIXED from Majid 2014.** van Hasselt 2013
  cites the upstream eribulin popPK model as "unpublished data on file, Dr Z.
  Hussein" (Eisai Limited, Hatfield, UK) and does not report PK parameter
  values in the paper. The Majid 2014 publication (J Clin Pharmacol
  54:1134-1143, doi:10.1002/jcph.315) is the published form of the same
  Eisai eribulin popPK analysis and shares the same first-author affiliation
  and Hussein co-authorship. The Majid 2014 PK parameters are used here
  verbatim from the encoding already in the registry at
  `Kawamura_2018_eribulin.R`. Operator decision recorded in sidecar
  request-001 of task
  frompeople-910-van_hasselt_2013_british_journal_of_clinical_ph.
- **PK dose-on-Q covariate NOT encoded.** van Hasselt 2013 mentions "the
  covariate dose was related to intercompartmental clearance" in the PK
  description, but the Majid 2014 published equations (reproduced in
  Kawamura 2018 section 2.3) do NOT include a dose-on-Q effect. Two
  possibilities: (a) the pre-2013 snapshot of the Eisai analysis referenced
  by van Hasselt 2013 had a dose-on-Q covariate that was dropped between
  that snapshot and the Majid 2014 publication, or (b) Kawamura 2018's
  reproduction of Majid 2014 is incomplete on this one covariate. The
  numerical dose-on-Q coefficient is not on disk in either paper. This
  model omits it; in practice, dose-on-Q has small impact on AUC for a
  fixed-dose regimen and does not affect the PD structure that is the
  primary focus of van Hasselt 2013.
- **CONMED_GCSF is a new canonical.** Following the Phase 3 covariate-column
  gate, `CONMED_GCSF` was registered as a new canonical binary covariate in
  `inst/references/covariate-columns.md` (operator decision recorded in
  sidecar request-002). Follows the well-established `CONMED_<drug>` pattern.
- **Table 5 incidence comparison caveat -- reference-covariate simulation
  under-represents the deep-nadir tail.** The paper's Table 5 (26.6% grade
  3, 16.85% grade 4) is derived from the no-covariate BASE model with
  n = 2000 stochastic simulations, where PD variability is absorbed into
  the IIV alone (base-model IIV on MTT = 23.2 CV%, on gamma = 19.9 CV%). The
  vignette simulates with the FINAL COVARIATE model at reference covariate
  values (n = 200), where the corresponding IIVs have been shrunk to 13.9
  and 12.2 CV% because part of the between-subject variability was absorbed
  into the retained ALB / ALP / TBILI / LDH / AST covariates. At reference
  covariates all covariate multipliers evaluate to 1, so the median simulated
  ANC trajectory matches the base-model prediction, but the tails are
  narrower. The vignette's simulated grade 3 = 23.5% is within one sample SE
  of the Table 5 base-model 26.6% (SE_n=200 ~ 3.1%). The vignette's
  simulated grade 4 = 4.0% is markedly lower than the Table 5 base-model
  16.85% (SE_n=200 ~ 2.6%); this ~13-percentage-point gap reflects the
  reference-covariate-averaging effect described in the paper's Discussion:
  the deepest grade 4 nadirs come from low-albumin patients (Table 5 albumin
  50% reduced row shows a +18.95 pp change on grade 4 for the SLOPE-albumin
  effect alone) and simulating the full cohort covariate distribution rather
  than fixing at the reference value would move the grade 4 incidence closer
  to Table 5. For a demonstration of that effect, users can substitute a
  covariate cohort drawn from Table 2's IQR distributions rather than the
  fixed-reference cohort used here.
- **No PKNCA comparison against published values.** van Hasselt 2013 does
  not tabulate NCA parameters for eribulin plasma concentrations. The PKNCA
  block is included as a sanity check on the Majid 2014 PK layer only, not
  as a comparison against the paper's reported values. For a comparison
  against Majid 2014's own NCA, see the future extraction task for that
  paper.
- **BSA distribution parameterisation.** The Methods section reports a
  simulation BSA distribution of mean = 1.57 m^2, sd = 0.22, truncated to
  [1, 3] m^2. Table 2 reports the observed cohort BSA median 1.73 m^2 (IQR
  1.61-1.86), a slightly right-skewed distribution. The vignette follows
  the Methods parameterisation because that is what van Hasselt 2013 used
  for its own simulation study.
- **Cohort size (n = 200) is smaller than the paper's n = 2000.** Reduced
  from the paper's simulation cohort size to keep the pkgdown vignette
  render under the 5-minute wall-clock gate. Monte-Carlo error on the
  grade-3 / grade-4 incidence estimates is proportionally larger but the
  point estimates remain within the expected range.
- **G-CSF prevalence 22.7%** matches Table 2 (382 / 1684 subjects). The
  paper does not describe temporal G-CSF administration patterns; the
  vignette assigns CONMED_GCSF as a per-subject baseline indicator, matching
  the paper's covariate encoding.
