---
title: "Gemcitabine + dFdU (Jiang 2008)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Gemcitabine + dFdU (Jiang 2008)}
  %\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: Jiang X, Galettis P, Links M, Mitchell PL, McLachlan AJ.
  Population pharmacokinetics of gemcitabine and its metabolite in
  patients with cancer: effect of oxaliplatin and infusion rate.
  *Br J Clin Pharmacol.* 2008;65(3):326-333.
* Article: <https://doi.org/10.1111/j.1365-2125.2007.03040.x>

The model is a joint two-compartment parent + two-compartment metabolite
population PK model fitted simultaneously to gemcitabine and its
inactive metabolite 2',2'-difluorodeoxyuridine (dFdU) in 94 adult
patients with cancer pooled from three Australian clinical trials
(Methods, page 327; Table 1). Gemcitabine is administered as a 30-min or
100-min IV infusion; the dose elimination from the gemcitabine central
compartment is treated as the apparent formation flux into the dFdU
central compartment because the fraction-of-conversion F is not
identifiable from plasma data (page 329, "The fraction of conversion
(F) from gemcitabine to dFdU could not be independently estimated and
hence apparent values of Q_dFdU/F, CL_dFdU/F, V_C,dFdU/F and V_P,dFdU/F
are reported"). Page 330 gives the final covariate equations:

```
CL_dFdU/F (L/min) = 0.04 * (1 + 0.48 * CRCL/70)
V_C,dFdU/F (L)    = 46 * (BSA/1.73)^0.93 * 0.65^GEMOX * 0.54^OXGEM * 1.24^NSCLC
```

# Population

The pooled cohort comprises 94 adult patients with cancer across three
clinical studies (Study 1: open-label Phase I dose-escalation, n=21,
diverse tumour types, gemcitabine-then-oxaliplatin; Study 2: randomised
Phase II sequencing trial, n=43, advanced NSCLC, gem-then-ox vs
ox-then-gem; Study 3: randomised crossover, n=31, diverse tumour types,
30-min vs 100-min gemcitabine infusion without oxaliplatin). Baseline
demographics from Table 2: age 33-85 years (median 65), body weight
37-120 kg (median 73), body surface area 1.2-2.5 m^2 (median 1.8),
estimated creatinine clearance 37-150 mL/min (median 83), 68 male and
26 female. Tumour-type split: NSCLC 47, pancreatic 14, diverse 33.
Oxaliplatin sequence split: 31 gemcitabine alone, 38 gemcitabine-
then-oxaliplatin, 25 oxaliplatin-then-gemcitabine. Concentration data:
652 plasma gemcitabine and 1130 plasma dFdU samples across 122 PK
profiles (HPLC-UV at 272 nm; LOQ 0.2 mg/L for both analytes; page 328).

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

# Source trace

The per-parameter origin is recorded as an in-file comment next to each
`ini()` entry in `inst/modeldb/specificDrugs/Jiang_2008_gemcitabine.R`.
The table below collects them in one place for review.

| Equation / parameter           | Value                       | Source location                                                            |
|--------------------------------|-----------------------------|----------------------------------------------------------------------------|
| `lcl`                          | log(2.7 L/min)              | Table 3 row CL_gem (estimate 2.7)                                          |
| `lvc`                          | log(15 L)                   | Table 3 row V_C,gem (estimate 15)                                          |
| `lq`                           | log(0.7 L/min)              | Table 3 row Q_gem (estimate 0.7)                                           |
| `lvp`                          | log(15 L)                   | Table 3 row V_P,gem (estimate 15)                                          |
| `lcl_dfdu`                     | log(0.04 L/min)             | Table 3 row CL_dFdU/F (estimate 0.04); page 330 covariate equation         |
| `lvc_dfdu`                     | log(46 L)                   | Table 3 row V_C,dFdU/F (estimate 46); page 330 covariate equation          |
| `lq_dfdu`                      | log(0.2 L/min)              | Table 3 row Q_dFdU/F (estimate 0.2)                                        |
| `lvp_dfdu`                     | log(192 L)                  | Table 3 row V_P,dFdU/F (estimate 192)                                      |
| `e_crcl_cl_dfdu`               | +0.48                       | Page 330 equation: CL_dFdU/F = 0.04 * (1 + 0.48 * CRCL/70)                 |
| `e_bsa_vc_dfdu`                | +0.93                       | Page 330 equation: V_C,dFdU/F = 46 * (BSA/1.73)^0.93 * ...                 |
| `e_gemox_vc_dfdu`              | 0.65                        | Page 330 equation: ... * 0.65^GEMOX                                        |
| `e_oxgem_vc_dfdu`              | 0.54                        | Page 330 equation: ... * 0.54^OXGEM                                        |
| `e_tumtp_nsclc_vc_dfdu`        | 1.24                        | Page 330 equation: ... * 1.24^NSCLC                                        |
| `etalcl`                       | omega^2 = 0.0918            | Table 3 BSV CL_gem 31% CV -> log(1 + 0.31^2)                               |
| `etalvc`                       | omega^2 = 0.1417            | Table 3 BSV V_C,gem 39% CV -> log(1 + 0.39^2)                              |
| `etalq`                        | omega^2 = 0.1770            | Table 3 BSV Q_gem 44% CV -> log(1 + 0.44^2)                                |
| `etalvp`                       | omega^2 = 0.3434            | Table 3 BSV V_P,gem 64% CV -> log(1 + 0.64^2)                              |
| `etalcl_dfdu`                  | omega^2 = 0.1155            | Table 3 BSV CL_dFdU/F 35% CV -> log(1 + 0.35^2)                            |
| `etalvc_dfdu`                  | omega^2 = 0.0223            | Table 3 BSV V_C,dFdU/F 15% CV -> log(1 + 0.15^2)                           |
| `etalq_dfdu`                   | omega^2 = 0.0808            | Table 3 BSV Q_dFdU/F 29% CV -> log(1 + 0.29^2)                             |
| `etalvp_dfdu`                  | omega^2 = 0.1349            | Table 3 BSV V_P,dFdU/F 38% CV -> log(1 + 0.38^2)                           |
| `propSd`                       | 0.40                        | Table 3 row "Residual variability sigma_GEM" (40% CV)                      |
| `propSd_dfdu`                  | 0.19                        | Table 3 row "Residual variability sigma_dFdU" (19% CV)                     |
| `d/dt(central)`                | -kel_gem * central - k12 * central + k21 * peripheral1 | Methods page 328 ("two-compartment ... FO elimination") + ADVAN6 parent-metabolite parameterisation |
| `d/dt(peripheral1)`            | k12 * central - k21 * peripheral1                      | Methods page 328                                                           |
| `d/dt(central_dfdu)`           | kel_gem * central - kel_dfdu * central_dfdu - k_dfdu_12 * central_dfdu + k_dfdu_21 * peripheral1_dfdu | Methods page 329 ("two-compartment ... FO formation rate") |
| `d/dt(peripheral1_dfdu)`       | k_dfdu_12 * central_dfdu - k_dfdu_21 * peripheral1_dfdu | Methods page 329                                                           |

# Virtual cohort

Original observed data are not publicly available. The cohort below
uses 300 virtual patients distributed across the three regimen groups
(100 gemcitabine alone, 100 gem-then-ox, 100 ox-then-gem) with body
surface area, creatinine clearance, and NSCLC indicator distributions
that bracket the published Table 2 demographics. Each patient receives
a single gemcitabine 1250 mg/m^2 30-min IV infusion (the modal Study 1
/ Study 2 dose); the absolute mg dose is computed per subject as
`1250 * BSA`.

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

n_per_grp <- 100L

make_cohort <- function(n, label, gemox, oxgem, p_nsclc, id_offset = 0L) {
  tibble::tibble(
    id          = id_offset + seq_len(n),
    BSA         = pmin(pmax(round(rlnorm(n, log(1.80), 0.10), 2), 1.2), 2.5),
    CRCL        = pmin(pmax(round(rlnorm(n, log(83),  0.30), 0), 37), 150),
    TUMTP_NSCLC = as.integer(runif(n) < p_nsclc),
    GEMOX       = gemox,
    OXGEM       = oxgem,
    treatment   = label,
    dose_mg_m2  = 1250
  )
}

# NSCLC prevalences by regimen track Table 1: Study 2 (gem-then-ox + ox-then-gem)
# was NSCLC-only (47/63 patients = 75% of combination arm); Study 1 (gem-then-ox
# in diverse tumour types) contributes 21 non-NSCLC patients to gem-then-ox.
# Study 3 (gemcitabine alone in diverse tumour types) is overwhelmingly
# non-NSCLC. Approximate prevalences below.
subj <- dplyr::bind_rows(
  make_cohort(n_per_grp, "gem alone",   gemox = 0L, oxgem = 0L, p_nsclc = 0.10, id_offset =   0L),
  make_cohort(n_per_grp, "gem-then-ox", gemox = 1L, oxgem = 0L, p_nsclc = 0.55, id_offset = 100L),
  make_cohort(n_per_grp, "ox-then-gem", gemox = 0L, oxgem = 1L, p_nsclc = 1.00, id_offset = 200L)
)

# Observation grid: dense early (capture gem infusion + post-infusion
# fall), then sparser out to ~1500 min so the dFdU terminal phase is
# resolved (matches Figure 6 time axis 0-2000 min).
obs_times <- c(
  seq(0,   60,   by =   2),   # infusion + first hour
  seq(75,  240,  by =  15),   # 1-4 h
  seq(270, 480,  by =  30),   # 4-8 h
  seq(540, 1500, by =  60)    # 9-25 h
)

build_events <- function(row) {
  amt_mg <- row$dose_mg_m2 * row$BSA
  dose   <- tibble::tibble(
    id = row$id, time = 0,
    amt = amt_mg, evid = 1L, cmt = "central", dur = 30
  )
  # Two observation streams: cmt = "Cc" (gemcitabine) and cmt = "Cc_dfdu"
  # (dFdU). rxSolve returns all model outputs at every observation row,
  # so a single stream is sufficient -- but we keep a single set of
  # observation rows to avoid duplicate (id, time, evid) keys.
  obs <- tibble::tibble(
    id = row$id, time = obs_times,
    amt = 0, evid = 0L, cmt = "Cc", dur = 0
  )
  dplyr::bind_rows(dose, obs) |>
    dplyr::mutate(
      BSA         = row$BSA,
      CRCL        = row$CRCL,
      TUMTP_NSCLC = row$TUMTP_NSCLC,
      GEMOX       = row$GEMOX,
      OXGEM       = row$OXGEM,
      treatment   = row$treatment
    )
}

events <- do.call(
  rbind,
  lapply(seq_len(nrow(subj)), function(i) build_events(subj[i, ]))
)
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
```

# Simulation

```{r simulate}
mod <- readModelDb("Jiang_2008_gemcitabine")

sim <- rxode2::rxSolve(
  mod,
  events = events,
  keep   = c("treatment", "BSA", "CRCL", "TUMTP_NSCLC", "GEMOX", "OXGEM")
) |>
  as.data.frame()
```

For deterministic replication of typical-value profiles (Figures 1-3
of the paper depict population predictions from the typical-value
parameters), zero out the random effects and re-solve one
representative subject per regimen.

```{r simulate-typical}
mod_typ <- mod |> rxode2::zeroRe()

typ_subj <- tibble::tribble(
  ~id, ~BSA, ~CRCL, ~TUMTP_NSCLC, ~GEMOX, ~OXGEM, ~treatment,
   1L,  1.8,    83,           0L,     0L,     0L, "gem alone",
   2L,  1.8,    83,           1L,     1L,     0L, "gem-then-ox",
   3L,  1.8,    83,           1L,     0L,     1L, "ox-then-gem"
)

events_typ <- do.call(
  rbind,
  lapply(seq_len(nrow(typ_subj)), function(i)
    build_events(typ_subj[i, ] |> dplyr::mutate(dose_mg_m2 = 1250))
  )
)

sim_typ <- rxode2::rxSolve(
  mod_typ,
  events = events_typ,
  keep   = c("treatment", "BSA", "CRCL", "TUMTP_NSCLC", "GEMOX", "OXGEM")
) |>
  as.data.frame()
```

# Replicate published figures

## Figure 5 -- gemcitabine concentration-time profile

Replicates Figure 5 of Jiang 2008 (predictive check of gemcitabine over
the 0-200 min window). The plot shows 5th, 50th, and 95th percentiles
of the stochastic simulation against time, by regimen. Gemcitabine
disposition is unaffected by the dFdU covariates (oxaliplatin order
and tumour-type indicators only affect V_C,dFdU/F), so all three
regimens should overlay closely.

```{r figure-5}
fig5 <- sim |>
  dplyr::filter(time > 0, time <= 200) |>
  dplyr::group_by(treatment, 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(fig5, aes(time, Q50, colour = treatment, fill = treatment)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.20, colour = NA) +
  geom_line(linewidth = 0.8) +
  scale_y_log10() +
  labs(
    x       = "Time (min)",
    y       = "Gemcitabine (mg/L)",
    title   = "Figure 5 -- gemcitabine predictive check",
    caption = "Replicates Figure 5 of Jiang 2008. Lines: 50th percentile; bands: 5th-95th percentiles."
  )
```

## Figure 6 -- dFdU concentration-time profile

Replicates Figure 6 of Jiang 2008 (dFdU predictive check, 0-2000 min).
The figure axes in the paper span up to 100 mg/L on a log scale and out
to ~1500 min in time. Unlike gemcitabine, dFdU disposition is sensitive
to the oxaliplatin-sequence indicators (V_C,dFdU/F decreased by 35%
when GEMOX=1, 46% when OXGEM=1) and the NSCLC indicator
(V_C,dFdU/F increased by 24% when NSCLC=1), so per-regimen overlays
should diverge in early-phase peak height.

```{r figure-6}
fig6 <- sim |>
  dplyr::filter(time > 0, time <= 1500) |>
  dplyr::group_by(treatment, time) |>
  dplyr::summarise(
    Q05 = quantile(Cc_dfdu, 0.05, na.rm = TRUE),
    Q50 = quantile(Cc_dfdu, 0.50, na.rm = TRUE),
    Q95 = quantile(Cc_dfdu, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(fig6, aes(time, Q50, colour = treatment, fill = treatment)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.20, colour = NA) +
  geom_line(linewidth = 0.8) +
  scale_y_log10() +
  labs(
    x       = "Time (min)",
    y       = "dFdU (mg/L)",
    title   = "Figure 6 -- dFdU predictive check",
    caption = "Replicates Figure 6 of Jiang 2008. Lines: 50th percentile; bands: 5th-95th percentiles."
  )
```

## Typical-value gemcitabine and dFdU profiles by regimen

The deterministic (zero-random-effects) simulation lets the covariate
effects on dFdU central volume be inspected without the noise from
between-subject variability.

```{r figure-typical}
typ_long <- sim_typ |>
  dplyr::filter(time > 0) |>
  dplyr::select(id, time, treatment, Cc, Cc_dfdu) |>
  tidyr::pivot_longer(
    cols      = c(Cc, Cc_dfdu),
    names_to  = "analyte",
    values_to = "conc"
  ) |>
  dplyr::mutate(
    analyte = factor(analyte, levels = c("Cc", "Cc_dfdu"),
                     labels = c("gemcitabine", "dFdU"))
  )

ggplot(typ_long, aes(time, conc, colour = treatment)) +
  geom_line(linewidth = 0.7) +
  facet_wrap(~analyte, scales = "free_y") +
  scale_y_log10() +
  labs(
    x       = "Time (min)",
    y       = "Concentration (mg/L)",
    title   = "Typical-value profiles by regimen",
    caption = "Zero-random-effects simulation: BSA = 1.8 m^2, CRCL = 83 mL/min, NSCLC = 0/1/1 across panels."
  )
```

# PKNCA validation

Two PKNCA passes are run -- one for gemcitabine `Cc` and one for the
dFdU metabolite `Cc_dfdu`. Each carries the `treatment` grouping so the
NCA summaries can be compared per regimen.

```{r pknca-gem}
gem_conc <- sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, treatment)

gem_conc <- dplyr::bind_rows(
  gem_conc,
  gem_conc |> 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 <- events |>
  dplyr::filter(evid == 1) |>
  dplyr::select(id, time, amt, treatment)

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

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

gem_nca <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(gem_conc_obj, dose_obj, intervals = intervals_gem)
)
```

```{r pknca-dfdu}
dfdu_conc <- sim |>
  dplyr::filter(!is.na(Cc_dfdu)) |>
  dplyr::select(id, time, Cc_dfdu, treatment) |>
  dplyr::rename(Cc = Cc_dfdu)

dfdu_conc <- dplyr::bind_rows(
  dfdu_conc,
  dfdu_conc |> dplyr::distinct(id, treatment) |>
    dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, treatment, time, .keep_all = TRUE) |>
  dplyr::arrange(id, treatment, time)

dfdu_conc_obj <- PKNCA::PKNCAconc(
  dfdu_conc, Cc ~ time | treatment + id,
  concu = "mg/L", timeu = "min"
)

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

dfdu_nca <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(dfdu_conc_obj, dose_obj, intervals = intervals_dfdu)
)
```

## Comparison against published parameter values

Jiang 2008 does not tabulate NCA summaries (Cmax, AUC, half-life) per
regimen -- the published comparison is instead a visual predictive
check (Figures 5 and 6) plus the per-parameter table (Table 3).
Below is a sanity check that the typical-value gemcitabine clearance
recovered by the simulation matches the published CL_gem = 2.7 L/min
within the expected NCA approximation uncertainty.

```{r sanity}
gem_res <- as.data.frame(gem_nca$result) |>
  dplyr::filter(PPTESTCD == "aucinf.obs") |>
  dplyr::mutate(dose = 1250 * 1.8) |>
  dplyr::group_by(treatment) |>
  dplyr::summarise(
    median_CL_app = median(dose / PPORRES, na.rm = TRUE),
    .groups       = "drop"
  )

knitr::kable(
  gem_res,
  caption = "Apparent gemcitabine clearance (= dose / AUCinf.obs) by regimen. Expected: CL_gem ~ 2.7 L/min (Jiang 2008 Table 3)."
)
```

# Assumptions and deviations

- **Between-occasion variability (BOV) is not encoded.** Jiang 2008
  Table 3 reports BOV on CL_gem (12% CV), V_C,gem (93% CV),
  CL_dFdU/F (11% CV), and V_C,dFdU/F (10% CV) in addition to the BSV
  values. nlmixr2 has no native occasion-level random-effect slot for
  a static model definition, so only BSV is encoded in `ini()`. Users
  fitting the model to multi-cycle data should consider adding an
  occasion-level eta layer externally.
- **NSCLC prevalence by regimen approximated from Table 1.** The
  paper does not publish a 3 x 2 contingency of regimen x NSCLC; the
  virtual cohort above uses Study 2 (NSCLC-only) + Study 1 (diverse
  tumour types, gem-then-ox) + Study 3 (diverse tumour types,
  gemcitabine alone) to set per-regimen NSCLC prevalences at 10% /
  55% / 100% respectively. The structural model is unaffected by
  this choice.
- **Demographics distributions approximated.** Body surface area and
  creatinine clearance distributions are sampled from log-normal
  curves whose median bracket Table 2's published medians; the
  actual cohort distribution is unknown from the paper.
- **Cmax / AUC NCA comparison.** Jiang 2008 reports a visual
  predictive check (Figures 5 and 6) and Table 3 parameter
  estimates rather than tabulated NCA values per regimen. The
  comparison section above therefore checks the apparent
  clearance recovered from the simulation against the published
  CL_gem = 2.7 L/min as a sanity check rather than rendering a
  full ncaComparisonTable.
- **dFdU `kel` is the apparent elimination rate.** The dFdU
  parameters (CL_dFdU/F, V_C,dFdU/F, etc.) are all apparent values
  because the fraction-of-conversion F from gemcitabine to dFdU is
  not identifiable (page 329). The simulated dFdU concentrations
  are nevertheless on the real plasma scale because the F factor
  cancels in `concentration = mass/F / volume/F = mass/volume`.
- **Allometric / size scaling not included.** Body weight was
  screened during covariate development but not retained (page 330,
  "No other covariates ... investigated in this study were
  identified as having a significant impact"). The final model
  uses BSA (on V_C,dFdU/F only) rather than weight as the body-size
  descriptor.
