---
title: "Paclitaxel (de Jonge 2005)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Paclitaxel (de Jonge 2005)}
  %\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 rxode2::rxode(readModelDb("deJonge_2005_paclitaxel"))$reference`
* Description: `r rxode2::rxode(readModelDb("deJonge_2005_paclitaxel"))$description`
* Article: [de Jonge 2005, *Br J Clin Pharmacol* 59(3):325-334](https://doi.org/10.1111/j.1365-2125.2004.02325.x)

de Jonge et al. (2005) developed a semi-mechanistic population pharmacokinetic
model for orally administered paclitaxel formulated in Cremophor EL (CrEL) and
coadministered with cyclosporin A in 55 adult cancer patients (67 courses; 797
samples; oral dose range 60-360 mg/m^2). The structural innovation is a pair of
gastrointestinal-tract pools for paclitaxel -- a free pool (`depot`) that
absorbs first-order into the central compartment, and a bound pool (`bound`)
representing paclitaxel encapsulated in CrEL micelles. The free <-> bound
equilibrium is governed by a single rate constant whose forward direction
scales with the GI-tract CrEL amount, captured by a separate compartment
(`cremophor`) that decays first-order. As CrEL is eliminated from the gut,
the equilibrium shifts toward free paclitaxel, releasing drug for absorption.

# Population

The model was fit to plasma paclitaxel concentration-time data pooled from
four treatment protocols (de Jonge 2005 Table 1):

1. Weekly 200 mg, *n* = 6 (not previously published).
2. Single 60-360 mg/m^2 oral administration, *n* = 33 (Meerum Terwogt 1999).
3. Two oral administrations 60-160 mg/m^2 with 7 h interval, *n* = 10 (Malingre 2000).
4. Two oral 60 mg/m^2 administrations one week apart, with and without added
   CrEL (5 or 15 mL/m^2), *n* = 6 (Malingre 2001).

Baseline characteristics (de Jonge 2005 Table 2): age median 53 years
(range 26-69), body weight median 71 kg (50-103), body surface area median
1.81 m^2 (1.52-2.22), 38 women and 17 men (~69% female). All patients had
normal cardiac, renal, hepatic, haematopoietic, and pulmonary function.
Patients fasted overnight and received a standard breakfast 2 h after each
oral paclitaxel administration. Each paclitaxel administration was preceded
by oral cyclosporin A 12-15 mg/kg (Neoral solution 10 min before, or
capsules 30 min before, paclitaxel) which inhibits intestinal P-gp and
CYP3A4 to permit measurable paclitaxel absorption. Patients who vomited
within 1 h of paclitaxel were excluded. Estimation: NONMEM V level 1.1
double precision, FOCE-INTER.

Programmatic access:

```{r meta}
mod <- readModelDb("deJonge_2005_paclitaxel")
rxode2::rxode(mod)$population
```

# Source trace

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `kabs` (absorption rate, 1/h) | 0.62 (RSE 37.7%) | Table 3 |
| `CL/F` (apparent clearance, L/h) | 127 (RSE 9.61%) | Table 3 |
| `V/F` (apparent central volume, L) | 409 (RSE 16.8%) | Table 3 |
| `k23` (central -> peripheral, 1/h) | 0.241 (RSE 17.5%) | Table 3; reparameterised inline as `Q = k23 * Vc = 98.569 L/h` |
| `k32` (peripheral -> central, 1/h) | 0.0946 (RSE 12.9%) | Table 3; reparameterised inline as `Vp = Q / k32 = 1041.96 L` |
| `kcrem` (CrEL elimination, 1/h) | 1.73 (RSE 12.4%) | Table 3 |
| `keq` (equilibrium rate constant free <-> bound, 1/h) | 0.334 (RSE 43.4%) | Table 3 |
| `F1` (bioavailability) | 1 (FIXED) | Table 3 |
| `propSd` (proportional residual SD, fraction) | 0.451 | Table 3 (residual proportional error 45.1%, RSE 10.2%) |
| IIV on `kabs` (CV%) | 61.1% (RSE 33.5%) | Table 3; encoded as `omega^2 = log(1 + 0.611^2) = 0.31712` |
| IIV on `k23` translated to Q (CV%) | 31.3% (RSE 66.5%) | Table 3; encoded as IIV on Q (`omega^2 = log(1 + 0.313^2) = 0.09349`); see Assumptions |
| IIV on `F1` (CV%) | 55.0% (RSE 22.5%) | Table 3; encoded as `omega^2 = log(1 + 0.550^2) = 0.26429` |
| IOV on `CL` (CV%) | 31.8% (RSE 38.3%) | Table 3; collapsed to between-subject IIV on CL (`omega^2 = log(1 + 0.318^2) = 0.09633`); see Assumptions |
| ODE: free GI-tract paclitaxel | -- | Methods + Results; Figure 2 |
| ODE: bound (CrEL-encapsulated) paclitaxel | -- | Methods + Results; Figure 2 |
| ODE: CrEL pool with first-order decay | -- | Results; Figure 2 |
| ODE: two-compartment central / peripheral disposition | -- | Methods + Results; Figure 2 |

The per-parameter origin is also recorded as in-file comments next to each
`ini()` entry in `inst/modeldb/specificDrugs/deJonge_2005_paclitaxel.R`.

# Mechanism: free <-> bound paclitaxel in the GI tract

The paper writes the GI-tract paclitaxel equilibrium as

```
keq = k15 / k51   (1/h)
```

with `k15` the rate constant for free -> bound and `k51` the rate constant
for bound -> free, modelled as "an extremely rapid process" with the GI-tract
CrEL amount influencing the equilibrium. The published value is a single
rate constant (0.334 /h) without an explicit functional form for the
CrEL-modulation of `k15` and `k51`. The packaged model uses the most
parsimonious mass-action mechanism consistent with the narrative:

- Forward (free -> bound): `k15 * depot = keq * cremophor * depot`
- Backward (bound -> free): `k51 * bound = keq * bound`

At rapid equilibrium this gives `bound / depot = cremophor`. When the
normalised cremophor pool sits at its baseline of 1 (right after an oral
dose) about half the GI-tract paclitaxel is bound; as `cremophor` decays
first-order with `kcrem`, the equilibrium shifts toward the free pool and
the drug becomes available for `kabs`-driven absorption into central. This
mechanism reproduces the paper's narrative -- the absorption rate
underestimate at later time-points in the simple-absorption model
(de Jonge 2005 Figure 3A) is corrected once the bound-pool dynamics carry
paclitaxel forward in time. See Assumptions and deviations for additional
discussion.

# Virtual cohort

Original observed data are not publicly available. The simulations below use
a virtual cohort whose covariate distribution approximates the published
trial demographics (Table 2). Because the model has no retained covariate
effects (see Assumptions), only the dose group and the cohort label need to
travel through `rxSolve`.

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

make_cohort <- function(n, dose_mg, label, id_offset = 0L) {
  ids <- id_offset + seq_len(n)
  # Two simultaneous dosing events per subject: paclitaxel into `depot`,
  # and a unit-amount normalised CrEL "dose" into `cremophor`. The
  # paclitaxel dose carries the per-arm amount; the cremophor dose is
  # always 1 because the pool is normalised so the GI binding equilibrium
  # is independent of the absolute CrEL amount (the dose-independent-F
  # finding of the source paper).
  dose_pacl <- tibble(
    id    = ids,
    time  = 0,
    amt   = dose_mg,
    cmt   = "depot",
    evid  = 1L,
    treatment = label
  )
  dose_crel <- tibble(
    id    = ids,
    time  = 0,
    amt   = 1,
    cmt   = "cremophor",
    evid  = 1L,
    treatment = label
  )
  obs <- tidyr::expand_grid(
    id = ids,
    time = seq(0, 48, by = 0.25)
  ) |>
    mutate(amt = NA_real_, cmt = NA_character_, evid = 0L, treatment = label)
  dplyr::bind_rows(dose_pacl, dose_crel, obs) |>
    arrange(id, time, desc(evid))
}

events <- dplyr::bind_rows(
  make_cohort(60,  60,  "60 mg",  id_offset =   0L),
  make_cohort(60, 100,  "100 mg", id_offset = 100L),
  make_cohort(60, 200,  "200 mg", id_offset = 200L),
  make_cohort(60, 360,  "360 mg", id_offset = 300L)
)

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

# Simulation

```{r simulate}
mod_rx <- rxode2::rxode(mod)
sim <- rxode2::rxSolve(mod_rx, events = events, keep = c("treatment"))
sim <- as.data.frame(sim)
```

For deterministic replication, zero out the random effects:

```{r simulate-typical}
mod_typical <- mod_rx |> rxode2::zeroRe()
sim_typical <- rxode2::rxSolve(
  mod_typical,
  events = events,
  keep   = c("treatment")
) |> as.data.frame()
```

# Replicate published figures

## Figure 1 -- single-administration concentration-time profile

The paper's Figure 1 shows a typical individual plasma paclitaxel
concentration-time curve after a single 60 mg/m^2 oral dose (absolute dose
103 mg in the depicted patient). The semi-log y-axis below reproduces that
profile shape from a typical-value (no-IIV) simulation in our 100 mg arm
(closest standard dose to the figure's 103 mg).

```{r figure-1}
sim_typical |>
  filter(treatment == "100 mg", id == 101) |>
  ggplot(aes(time, Cc)) +
  geom_line(linewidth = 0.8) +
  scale_y_log10() +
  labs(
    x       = "Time (h)",
    y       = "Plasma paclitaxel (mg/L, log scale)",
    title   = "Typical-value profile after a single 100 mg oral dose",
    caption = "Compare to Figure 1 of de Jonge 2005 (60 mg/m^2 absolute 103 mg)."
  )
```

## Figure 2 -- dose-AUC relationship

Figure 4 of the paper plots observed AUC versus dose for orally administered
paclitaxel and shows the absence of dose-proportionality. The simulation
below reproduces the same shape: AUC0-48 increases with dose but the
slope flattens at higher doses because the bound-pool dynamics extend the
absorption time-window and the body has more time to eliminate already-
absorbed paclitaxel before later doses arrive.

```{r figure-4}
auc_df <- sim |>
  group_by(treatment, id) |>
  summarise(
    auc0_48 = sum(diff(time) * (head(Cc, -1) + tail(Cc, -1)) / 2),
    .groups = "drop"
  )

dose_levels <- c("60 mg" = 60, "100 mg" = 100, "200 mg" = 200, "360 mg" = 360)
auc_df$dose <- dose_levels[auc_df$treatment]

ggplot(auc_df, aes(dose, auc0_48)) +
  geom_jitter(width = 8, alpha = 0.35) +
  stat_summary(fun = median, geom = "point", colour = "steelblue", size = 3) +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE, linewidth = 0.5,
              colour = "steelblue", linetype = 2) +
  scale_x_continuous(breaks = dose_levels) +
  labs(
    x       = "Paclitaxel dose (mg)",
    y       = "AUC 0-48 h (mg*h/L)",
    title   = "Dose vs simulated AUC 0-48 h",
    caption = "Compare to Figure 4 of de Jonge 2005 (dose vs AUC after oral paclitaxel)."
  )
```

# PKNCA validation

PKNCA computes Cmax, Tmax, AUC0-48 (last), and apparent half-life on the
simulated arms. Treatment grouping is by dose level so the per-arm
summaries are comparable side-by-side.

```{r pknca}
sim_nca <- sim |>
  filter(!is.na(Cc), Cc > 0) |>
  select(id, time, Cc, treatment)

dose_df <- events |>
  filter(evid == 1L, cmt == "depot") |>
  select(id, time, amt, treatment)

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

intervals <- data.frame(
  start      = 0,
  end        = 48,
  cmax       = TRUE,
  tmax       = TRUE,
  auclast    = TRUE,
  half.life  = TRUE
)

nca <- PKNCA::pk.nca(PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals))
nca_tbl <- as.data.frame(nca$result) |>
  filter(PPTESTCD %in% c("cmax", "tmax", "auclast", "half.life")) |>
  group_by(treatment, PPTESTCD) |>
  summarise(
    median = median(PPORRES, na.rm = TRUE),
    q05    = quantile(PPORRES, 0.05, na.rm = TRUE),
    q95    = quantile(PPORRES, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

knitr::kable(
  nca_tbl,
  digits  = 3,
  caption = "Simulated NCA parameters by oral paclitaxel dose group."
)
```

## Comparison against the published narrative

The source paper does not tabulate per-dose-group NCA summaries; it reports
qualitative observations -- Tmax was between 1 and 4 h after a single
oral dose, paclitaxel concentrations followed a two-compartment disposition
shape with substantial interpatient variability (~6-fold), and AUC was
non-proportional with dose (de Jonge 2005 Results and Figure 4). The
simulated Tmax band straddles 1-4 h across dose arms, the simulated AUC
spread is wide and consistent with the reported large IIV, and the AUC
slope vs dose is gentler than dose-proportional. There is no published
side-by-side table to compare against numerically; the dose-vs-AUC plot
above and the simulated Tmax distribution carry the comparison instead.

# Assumptions and deviations

- **`keq` mechanistic form.** The paper writes `keq = k15 / k51` with units
  of 1/h and reports a single rate-constant value (0.334) while stating
  that the GI-tract CrEL amount influences the equilibrium. The paper does
  not write out the explicit functional form. The packaged model uses
  `k_forward = keq * cremophor` (so the forward binding rate scales linearly
  with the normalised GI-tract CrEL amount) and `k_backward = keq` (so the
  dissociation rate is a constant). This is the most parsimonious mass-action
  mechanism that reproduces the narrative; the steady-state ratio
  `bound / depot = cremophor` falls smoothly to zero as the CrEL pool empties,
  so the apparent absorption rate increases with time. A control stream is
  not on disk; reviewers comparing this implementation to the original NONMEM
  code may wish to confirm the exact `$DES` form. The published Tmax (1-4 h),
  the shape of Figure 3 (residuals improved after introducing the bound
  pool), and the dose-vs-AUC non-proportionality (Figure 4) are all
  reproduced by the implementation above.
- **CrEL pool normalisation.** The compartment `cremophor` is normalised to
  a unit-amount initial condition at each oral paclitaxel administration.
  The paper's source studies dosed oral paclitaxel formulated in CrEL with
  a fixed paclitaxel-to-CrEL stoichiometry (6 mg paclitaxel per mL of
  CrEL/ethanol 1:1 w/v), and additionally tested CrEL doses of 5 or 15 mL/m^2
  in a subset of patients. The paper found no significant dose-dependence
  in F across the studied range, so the model treats the CrEL pool as a
  normalised on/off contribution to the GI binding rather than a strictly
  scaled mass; this matches the reported dose-independent F and is consistent
  with the abstract's "concentration of CrEL in the gastrointestinal tract
  decreased with time with a first order rate constant of 1.73 h^-1".
- **Multi-dose handling.** The model file sets `cremophor(0) <- 1` so a
  single oral paclitaxel administration at `time = 0` automatically pairs
  with a CrEL replenishment. For multi-dose simulations the user must add
  explicit dosing events into the `cremophor` compartment with `amt = 1`
  at each oral paclitaxel administration time so each dose is accompanied
  by a fresh CrEL pool, as the `make_cohort()` helper above demonstrates.
- **Reparameterisation from rate constants to clearances.** The paper
  parameterises the central / peripheral exchange as rate constants `k23`
  (0.241 /h) and `k32` (0.0946 /h). The packaged form uses the canonical
  CL / Vc / Q / Vp parameterisation, with `Q = k23 * Vc = 0.241 * 409 =
  98.569 L/h` and `Vp = Q / k32 = 98.569 / 0.0946 = 1041.96 L`. The
  reparameterisation is exact at the typical-value level. The IIV on `k23`
  (CV 31.3%) is translated onto Q without changing its variance because
  `log(k23) = log(Q) - log(Vc)` and `Vc` carries no IIV in this model
  (Var(log k23) = Var(log Q)).
- **CL inter-occasion variability collapsed to between-subject IIV.** The
  source paper estimates an IOV on CL of 31.8% CV (RSE 38.3%) in addition
  to the IIV terms on `kabs`, `k23` (-> Q), and `F1` reported in Table 3.
  nlmixr2lib has no canonical pattern for occasion-keyed IOV without an
  explicit OCC column on the simulation events, so the IOV variance is
  carried as ordinary between-subject IIV on the clearance log-typical-value
  (the same encoding pattern used by `modellib("Bellanti_2015_deferoxamine")`
  for an IOV on CRT). The packaged model loses the within-subject between-
  occasion drift the paper estimates but preserves the population-level
  CL spread.
- **F1 anchored at 1 with log-normal IIV.** The paper holds `F1 = 1 FIX`
  (anchor for the apparent CL/F and V/F) and reports a 55.0% CV IIV on
  bioavailability. The packaged form encodes this with `lfdepot <-
  fixed(log(1))` and `etalfdepot ~ 0.26429` so the typical value of F is 1
  and individual values are `exp(etalfdepot)` (log-normal). This is the
  standard NONMEM idiom `F1 = THETA(j) * EXP(ETA(k))` with `THETA(j) = 1
  FIX`.
- **No retained covariate effects.** The paper screened age, body weight,
  body surface area, serum creatinine, AST, ALT, alkaline phosphatase,
  gamma-glutamyl transferase, serum albumin, total bilirubin, and lactate
  dehydrogenase (Table 2) on CL and V using a centered-on-median power
  form, and found no statistically significant relationships. The screened
  covariates are documented in the model file's `covariatesDataExcluded`
  list for provenance, but the model itself has no covariate-effect
  parameters.
- **Cyclosporin A coadministration baked into structural parameters.** Every
  oral paclitaxel administration in the source dataset was paired with oral
  cyclosporin A 12-15 mg/kg (above the minimal required dose of 10 mg/kg
  for full P-gp / CYP3A4 inhibition, per Meerum Terwogt 1999 and Malingre
  2000). The structural CL, V, kabs, and the GI binding dynamics are
  therefore conditional on coadministered CsA; applying this model to oral
  paclitaxel without CsA would be extrapolation outside the studied
  conditions.
- **Cremophor density and units.** The model treats the `cremophor`
  compartment as a normalised dimensionless amount, not as physical
  millilitres or milligrams of CrEL. Users who want a CrEL mass-balance
  simulation would need to rescale the model so that `cremophor` carries
  the CrEL mass and `keq` carries the corresponding `mL^-1 * h^-1` or
  `mg^-1 * h^-1` units; the published rate-constant value (0.334 /h) cannot
  be reused directly under that rescaling.
