---
title: "Topotecan (Leger 2004)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Topotecan (Leger 2004)}
  %\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("Leger_2004_topotecan"))$meta
```

* Citation: `r mod_meta$reference`
* Description: `r mod_meta$description`
* Article (DOI): <https://doi.org/10.1038/sj.bjc.6601469>

This vignette validates the packaged `Leger_2004_topotecan` model -- a
two-compartment population PK model for oral and intravenous topotecan in
190 adult cancer patients pooled across five clinical trials and a drug
monitoring cohort -- against the source publication's Table 1 (baseline
demographics), Table 3 (final-model parameter estimates and inter-individual
variability), and the Results / Discussion narrative on inter-occasion
variability, oral bioavailability, and the limited-sampling strategy for
oral topotecan AUC.

# Population

Leger 2004 pooled topotecan plasma concentration data from 190 adult cancer
patients enrolled in five separate clinical trials (173 patients) plus 17
additional patients receiving therapeutic drug monitoring for other reasons.
Patients received topotecan either as a 30-minute intravenous infusion
(n = 72) at 0.2-2.4 mg/m^2/day, or as oral gelatin capsules (n = 118) at
0.15-2.7 mg/m^2/day, for 5-21 consecutive days. The cohort is European
(Toulouse, Institut Claudius-Regaud and Rotterdam, Erasmus MC -- Daniel den
Hoed Cancer Center). Median age was 55 years (range 18-76), median body
weight 70 kg (range 42-117), and median body surface area 1.80 m^2 (range
1.36-2.44, DuBois formula). Mean Cockcroft-Gault creatinine clearance was
80 mL/min (range 33-167) and mean serum creatinine was 87 umol/L (range
41-162). WHO performance status (PS) was 0, 1, 2, or 3 in 79, 98, 11, and
2 patients respectively. 141 of 190 patients had received previous
chemotherapy, and cisplatin pre-treatment was a feature of 115 patients
(prior regimen) and 30 patients (on Day 1 of topotecan). The total number
of analyzable plasma samples was 2064 (2-21 samples per patient, median 15
in cycle 1).

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/Leger_2004_topotecan.R`. The
table below collects them in one place; all values come from Leger 2004
Table 3 (final-covariate-model, cycle 1) unless otherwise noted.

| Parameter / equation                              | Value             | Source location |
|---------------------------------------------------|-------------------|-----------------|
| `lcl` (CL intercept; non-renal CL)                | log(12.8)         | Table 3, theta1 = 12.8 (95% CI 4.8) |
| `e_crcl_cl` (slope of CL on CrCl in L/h)          | 2.1               | Table 3, theta2 = 2.1 (95% CI 1.0) |
| `e_who_ps_cl` (fractional CL reduction per WHO PS)| 0.12              | Table 3, theta3 = 0.12 (95% CI 0.09) |
| `lvc` (slope of V1 on body weight, L/kg)          | log(0.58)         | Table 3, theta4 = 0.58 (95% CI 0.13) |
| `lvp` (V2, L)                                     | log(45.5)         | Table 3, theta5 = 45.5 (95% CI 7.0) |
| `lq` (Q, L/h)                                     | log(49.2)         | Table 3, theta6 = 49.2 (95% CI 16.9) |
| `lfdepot` (oral bioavailability F)                | log(0.324)        | Table 3, theta7 = 32.4% (95% CI 3.9) |
| `lka` (first-order absorption rate, 1/h)          | log(1.7)          | Table 3, theta9 = 1.7 (95% CI 0.6) |
| `ltlag` (oral absorption lag time, h)             | log(0.17)         | Table 3, theta10 = 0.17 (95% CI 0.03) |
| `etalcl ~ 0.08075`                                | log(0.29^2 + 1)   | Table 3, CL %CV = 29 (interindividual) |
| `etalvc ~ 0.14157`                                | log(0.39^2 + 1)   | Table 3, V1 %CV = 39 |
| `etalvp ~ 0.19193`                                | log(0.46^2 + 1)   | Table 3, V2 %CV = 46 |
| `etalq  ~ 0.49470`                                | log(0.80^2 + 1)   | Table 3, Q  %CV = 80 |
| `etalfdepot ~ 0.04727`                            | log(0.22^2 + 1)   | Table 3, F  %CV = 22 |
| `etalka ~ 0.34335`                                | log(0.64^2 + 1)   | Table 3, Ka %CV = 64 |
| `etaltlag ~ 0.01676`                              | log(0.13^2 + 1)   | Table 3, Lag %CV = 13 |
| `propSd <- 0.17`                                  | 0.17              | Results, oral-route proportional residual = 17% |
| `addSd  <- 0.09`                                  | 0.09 ug/L         | Results, oral-route additive residual = 0.09 ug/L |
| `cl <- (exp(lcl) + e_crcl_cl * crcl_Lh) * (1 - e_who_ps_cl * WHO_PS) * exp(etalcl)` | n/a | Table 3 covariate model row "CL = (theta1 + theta2 * CrCl) * (1 - theta3 * PS)"; CrCl converted from mL/min to L/h via 0.06 in model() per Table 3 footnote |
| `vc <- exp(lvc + etalvc) * WT`                    | n/a               | Table 3 row "V1 = theta4 * body weight" |
| `f(depot) <- exp(lfdepot + etalfdepot)`           | n/a               | Results, oral bioavailability mean 32.4%, CV 22% |
| `alag(depot) <- exp(ltlag + etaltlag)`            | n/a               | Results / Table 3, lag time mean 0.17 h, CV 13% |
| `Cc ~ add(addSd) + prop(propSd)`                  | n/a               | Results, combined error model (additive + proportional); oral-route values used (see Assumptions) |

# Virtual cohort

Original observed concentrations are not publicly available. The virtual
cohorts below approximate the published demographics: body weight log-normal
around the 70-kg median (38-125 kg range); Cockcroft-Gault creatinine
clearance log-normal around the 80 mL/min mean (truncated to the 33-167
mL/min range reported in Table 1); WHO performance status drawn from the
empirical 79 / 98 / 11 / 2 distribution. Two parallel cohorts are
constructed -- 60 oral subjects (1.5 mg/m^2 single dose) and 60 IV subjects
(1.5 mg/m^2 over 30 min) -- so the oral-versus-IV bioavailability check
later in this vignette has comparable sample sizes.

```{r cohort}
set.seed(20260610)

draw_cohort <- function(n, id_offset = 0L) {
  wt   <- exp(rnorm(n, mean = log(70),  sd = log(125 / 38) / 4))
  wt   <- pmin(pmax(wt, 42), 117)
  bsa  <- 0.007184 * wt^0.425 * 170^0.725
  bsa  <- pmin(pmax(bsa, 1.36), 2.44)
  crcl <- exp(rnorm(n, mean = log(80),  sd = log(167 / 33) / 4))
  crcl <- pmin(pmax(crcl, 33), 167)
  ps   <- sample(0:3, n, replace = TRUE,
                 prob = c(79, 98, 11, 2) / 190)
  tibble::tibble(
    id     = id_offset + seq_len(n),
    WT     = wt,
    BSA    = bsa,
    CRCL   = crcl,
    WHO_PS = ps
  )
}

n_oral <- 60L
n_iv   <- 60L

cov_oral <- draw_cohort(n_oral, id_offset =  0L)  |> mutate(treatment = "oral")
cov_iv   <- draw_cohort(n_iv,   id_offset = 1000L) |> mutate(treatment = "iv")

# Dosing target: 1.5 mg/m^2 single dose, ~the i.v. recommended dose in von
# Pawel 2001 cited in the Discussion. Oral subjects receive an oral capsule
# (depot compartment). IV subjects receive a 30-minute infusion (central
# compartment) per the protocol.
dose_mg <- function(bsa) round(1.5 * bsa, 3)

sample_times <- c(0, 0.25, 0.5, 0.75, 1, 1.5, 2, 3, 4, 6, 8, 12, 24)

make_subject <- function(row) {
  amt <- dose_mg(row$BSA)
  if (identical(row$treatment, "oral")) {
    dose_event <- tibble::tibble(
      id   = row$id,         time = 0,
      evid = 1L,             amt  = amt,
      rate = NA_real_,       cmt  = "depot"
    )
  } else {
    dose_event <- tibble::tibble(
      id   = row$id,         time = 0,
      evid = 1L,             amt  = amt,
      rate = amt / 0.5,      cmt  = "central"
    )
  }
  obs_events <- tibble::tibble(
    id   = row$id,           time = sample_times,
    evid = 0L,               amt  = NA_real_,
    rate = NA_real_,         cmt  = NA_character_
  )
  dplyr::bind_rows(dose_event, obs_events) |>
    dplyr::mutate(
      WT        = row$WT,
      CRCL      = row$CRCL,
      WHO_PS    = row$WHO_PS,
      treatment = row$treatment
    ) |>
    dplyr::arrange(time, dplyr::desc(evid))
}

build_events <- function(cov_tbl) {
  dplyr::bind_rows(lapply(seq_len(nrow(cov_tbl)), function(i) {
    make_subject(cov_tbl[i, ])
  }))
}

events_oral <- build_events(cov_oral)
events_iv   <- build_events(cov_iv)
events      <- dplyr::bind_rows(events_oral, events_iv)

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

# Simulation

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

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

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

# Replicate published figures

Leger 2004 has two main figures: Figure 1 shows two representative oral
subjects with limited (-9%) versus large (+104%) inter-day change in AUC;
Figure 2 shows the limited-sampling strategy ratio plot. Both depend on
individual observed data that is not in the public domain. The chunks below
instead reproduce the published model's typical-value oral and IV
concentration-time profiles, which faithfully exercise the structural
two-compartment model with first-order absorption and lag time.

## Typical-value oral and IV concentration-time profiles

```{r typical-profiles}
typ_summary <- sim_typical |>
  dplyr::filter(time > 0) |>
  dplyr::group_by(time, treatment) |>
  dplyr::summarise(
    Cc_median = median(Cc, na.rm = TRUE),
    .groups   = "drop"
  )

ggplot(typ_summary, aes(time, Cc_median, colour = treatment)) +
  geom_line(linewidth = 1) +
  scale_y_log10() +
  labs(
    x        = "Time (h)",
    y        = "Topotecan plasma concentration (ug/L, log scale)",
    colour   = "Route",
    title    = "Typical-value topotecan concentration-time profile after 1.5 mg/m^2",
    subtitle = "Two-compartment model with first-order absorption (oral) or 30-min infusion (IV)"
  ) +
  theme_minimal()
```

## Stochastic VPC over the oral cohort

```{r vpc-oral}
sim_stoch |>
  dplyr::filter(time > 0, treatment == "oral") |>
  dplyr::group_by(time) |>
  dplyr::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)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95),
              fill = "gray70", alpha = 0.5) +
  geom_line(linewidth = 0.9) +
  scale_y_log10() +
  labs(
    x        = "Time (h)",
    y        = "Topotecan plasma concentration (ug/L, log scale)",
    title    = "Stochastic VPC over the oral cohort (1.5 mg/m^2 single dose)",
    subtitle = "Median (line) and 90% prediction interval (ribbon)"
  ) +
  theme_minimal()
```

# PKNCA validation

Leger 2004 reports an observed mean bioavailability of 32.4% (95% CI +/- 3.9)
for oral topotecan and an inter-individual variability on dose-normalized
AUC of 4.8-fold for IV data and 7.6-fold for oral data (Discussion, p. 346).
The simulated NCA below compares oral-versus-IV AUC0-Inf to recover the
bioavailability and AUC variability.

```{r pknca}
sim_for_nca <- sim_stoch |>
  dplyr::filter(!is.na(Cc), time > 0) |>
  dplyr::select(id, time, Cc, treatment) |>
  as.data.frame()

doses_for_nca <- events |>
  dplyr::filter(evid == 1L) |>
  dplyr::select(id, time, amt, treatment) |>
  as.data.frame()

conc_obj <- PKNCA::PKNCAconc(
  data    = sim_for_nca,
  formula = Cc ~ time | treatment + id,
  concu   = "ug/L",
  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,
  half.life  = TRUE
)

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

nca_tab <- as.data.frame(nca_res$result) |>
  dplyr::select(treatment, id, start, end, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES)

knitr::kable(
  head(nca_tab, 12),
  digits  = 3,
  caption = "Simulated NCA per subject (first 12 rows)."
)
```

## Comparison against published bioavailability

```{r f-check}
auc_per_subject <- nca_tab |>
  dplyr::select(id, treatment, aucinf.obs) |>
  dplyr::left_join(
    dplyr::bind_rows(cov_oral, cov_iv) |>
      dplyr::select(id, BSA),
    by = "id"
  ) |>
  dplyr::mutate(
    dose_mg            = 1.5 * BSA,
    auc_per_mg_per_m2  = aucinf.obs / 1.5
  )

mean_auc_per_mg_m2 <- auc_per_subject |>
  dplyr::group_by(treatment) |>
  dplyr::summarise(
    mean_auc = mean(auc_per_mg_per_m2, na.rm = TRUE),
    sd_auc   = sd(auc_per_mg_per_m2,   na.rm = TRUE),
    auc_fold = max(auc_per_mg_per_m2, na.rm = TRUE) /
               min(auc_per_mg_per_m2, na.rm = TRUE),
    .groups  = "drop"
  )

f_observed <- mean_auc_per_mg_m2$mean_auc[mean_auc_per_mg_m2$treatment == "oral"] /
              mean_auc_per_mg_m2$mean_auc[mean_auc_per_mg_m2$treatment == "iv"]

knitr::kable(
  mean_auc_per_mg_m2,
  digits  = 2,
  caption = paste0("Dose-normalised AUC0-Inf by route. Empirical bioavailability ",
                   "(oral AUC / IV AUC) = ",
                   round(100 * f_observed, 1), "% (published value 32.4%).")
)
```

The empirical bioavailability recovered from the dose-normalised
oral-versus-IV AUC ratio should be close to the published 32.4%; small
discrepancies (a few percentage points) arise from the finite cohort size
and the log-normal IIV on F (22% CV) plus IIV on CL (29% CV) interacting
in the simulated NCA. The AUC fold-spread per route also varies with the
finite N. Leger 2004 reports 4.8-fold variability for IV and 7.6-fold for
oral.

# Assumptions and deviations

- **Residual error encoded for the oral route only.** Leger 2004 reports
  two route-specific residual error models in the Results section:
  proportional 11% with additive 0.64 ug/L for the IV data, and
  proportional 17% with additive 0.09 ug/L for the oral data. The
  different additive values largely reflect the different HPLC lower
  limits of quantification at the two analytical sites (Toulouse 0.5
  ng/mL for IV; Rotterdam 0.1 ng/mL for oral). The packaged model uses
  the oral-route residual error because (a) the paper's primary
  contribution is oral topotecan PK and the limited-sampling strategy is
  developed for oral, and (b) 118 of 190 patients (62%) are oral. For an
  IV-only simulation, override `propSd` to 0.11 and `addSd` to 0.64 ug/L
  in the loaded model before calling `rxSolve()`.
- **Interoccasion variability omitted.** Leger 2004 reports "interday
  variability" of 18% on CL, 49% on V1, and 28% on F as a separate IOV
  term in Table 3. The packaged model omits this layer because (a) the
  standard popPK extraction convention in this library is to capture only
  inter-individual variability and (b) IOV requires an OCC column in the
  event dataset that would constrain downstream use. The 18% / 49% / 28%
  IOV values are noted in the model file as comments so a user wishing
  to add an OCC dimension can wire them in.
- **No race / ethnicity covariate.** Leger 2004 does not report the
  cohort's race or ethnicity distribution (the paper notes only "French
  and Dutch oncology centres"); the virtual cohort therefore does not
  stratify on race. The `population$race_ethnicity` field records this
  as "Not reported".
- **Cohort body-surface area derived from weight only.** The virtual
  cohort uses the DuBois formula on body weight assuming a fixed 170 cm
  reference height to derive BSA for dosing (`1.5 mg/m^2 x BSA`). Real
  cohorts have independent height variation; this assumption is
  immaterial for the structural PK validation but means the simulated
  per-subject dose has a narrower distribution than the original cohort.
- **WHO PS draw uses the empirical 79 / 98 / 11 / 2 distribution.** The
  cohort's PS distribution is encoded directly (Table 1). PS values are
  treated as continuous-integer per the encoded covariate model
  (canonical column `WHO_PS`, added to `inst/references/covariate-columns.md`
  as part of this extraction). The `(1 - 0.12 * WHO_PS)` form means a
  WHO_PS = 3 patient has 36% lower CL than a reference WHO_PS = 0
  patient.
- **Cycle 1 vs Cycle 2 estimates.** The model file encodes the cycle-1
  final-model values from Table 3 (theta1 = 12.8, theta2 = 2.1, theta3 =
  0.12, theta4 = 0.58 etc.). Leger 2004 also fits a cycle-1-plus-cycle-2
  pooled model that gives slightly different coefficients (CL = (11.4 +
  2.5 * CrCl) * (1 - 0.06 * PS), V1 = 0.50 * body weight) plus
  inter-cycle variability of 18% on CL, 31% on V1, and 22% on F. These
  pooled-cycle values are not encoded; the paper highlights the cycle-1
  fit as the primary covariate model.
- **Alternative covariate models from Table 3 not encoded.** Table 3
  reports several reduced-form alternative covariate models (CL = theta1
  alone, CL = theta1 + theta2 * CrCl alone, CL = theta1 * (1 - theta3 *
  PS) alone, V1 = theta4 alone) with their objective-function deltas;
  these are model-development byproducts shown for comparison, not the
  final reported model, and are not encoded per the
  replicate-author-structure policy.
