---
title: "Ciclosporin (Wilhelm 2012)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Ciclosporin (Wilhelm 2012)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(nlmixr2lib)
library(PKNCA)
library(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
mod_meta <- rxode2::rxode(readModelDb("Wilhelm_2012_ciclosporin"))
```

# Model and source

* Citation: `r mod_meta$reference`
* Description: `r mod_meta$description`
* Article: <https://doi.org/10.1111/j.1365-2125.2011.04116.x>

# Population

Wilhelm 2012 enrolled 20 adults undergoing allogeneic haematopoietic stem
cell transplantation (HSCT) at VU University Medical Center, Amsterdam,
between January 2005 and February 2008. Demographics (Table 1): median
age 54 years (range 37-66), median body weight 84 kg (range 53-110),
median body surface area 2.02 m^2 (1.48-2.43), 13 male / 7 female. The
underlying haematological malignancies were acute myeloid leukaemia
(n=7), non-Hodgkin lymphoma (n=6), chronic lymphoblastic leukaemia
(n=2) and other diagnoses (n=5). Conditioning regimens were
fludarabine + cyclophosphamide (n=10), fludarabine + total-body
irradiation (n=6), or cyclophosphamide + total-body irradiation (n=4).
All subjects received fluconazole 50 mg once daily as routine
antimycotic prophylaxis throughout the sampling period; fluconazole is
a moderate CYP3A4 inhibitor that the authors highlight as the most
likely driver of the relatively low clearance reported here compared
with non-HSCT renal/liver transplant cohorts (Discussion, page
referencing Schultz et al. and Dotti et al.).

Ciclosporin was administered as a 2.5 mg/kg intravenous infusion over
3 h at the start of the conditioning scheme. After tolerance to oral
intake was established, dosing was switched to twice-daily oral
Neoral microemulsion at clinician-adjusted doses (mean 2.77 +/- 0.81
mg/kg per dose during the recorded oral profile) targeting trough
concentrations of 200-400 ug/L. Whole-blood ciclosporin was assayed by
fluorescence polarization immunoassay (FPIA, AxSYM Abbott) with LLOQ
80 ug/L and intra-assay CV below 10%; 436 concentration measurements
were available for the popPK analysis.

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

# Source trace

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

| Equation / parameter | Value | Source location (Wilhelm 2012) |
|----------------------|-------|--------------------------------|
| Two-compartment open model with first-order oral absorption + lag time | n/a | Methods 'Basic pharmacokinetic model'; Results Table 2 |
| IV dosing: 3 h infusion into central compartment | n/a | Methods 'Basic pharmacokinetic model' ('I.v. administration was modelled as a 3 h infusion in the central compartment') |
| Combined additive + proportional residual error | n/a | Methods 'Basic pharmacokinetic model' |
| `lcl` (CL) | log(21.9) L/h | Table 3, RSD 5.2% |
| `lvc` (V1) | log(16.6) L | Table 3, RSD 8.7% (abstract reports 18.3 L matching the bootstrap median ~18; the formal final estimate from Table 3 is 16.6 L) |
| `lq` (Q) | log(24.2) L/h | Table 3, RSD 9.3% |
| `lvp` (V2) | log(59.0) L | Table 3, RSD 8.8% |
| `lka` (ka) | log(0.280) 1/h | Table 3, RSD 14.6% |
| `lfdepot` (F) | log(0.710) | Table 3, RSD 9.9% |
| `ltlag` (tlag) | log(0.440) h | Table 3, RSD 5.5% |
| `etalcl` variance | log(1 + 0.222^2) = 0.04812 | Table 3 IIV CL 22.2% CV (RSD 55%) |
| `etalvc` variance | log(1 + 0.269^2) = 0.06988 | Table 3 IIV V1 26.9% CV (RSD 53%) |
| `etalq` variance | log(1 + 0.282^2) = 0.07650 | Table 3 IIV Q 28.2% CV (RSD 73%) |
| `etalvp` variance | log(1 + 0.306^2) = 0.08947 | Table 3 IIV V2 30.6% CV (RSD 62%) |
| `etalka` variance | log(1 + 0.438^2) = 0.17559 | Table 3 IIV ka 43.8% CV (RSD 66%) |
| `etalfdepot` variance | log(1 + 0.250^2) = 0.06062 | Table 3 IIV F 25.0% CV (RSD 64%) |
| `etaltlag` variance | log(1 + 0.181^2) = 0.03224 | Table 3 IIV tlag 18.1% CV (RSD 90%) |
| `addSd` (ug/L) | 65 | Table 3 'Additive error (ug/L)' = 65 (RSD 86%) |
| `propSd` (fraction) | 0.088 | Table 3 'Proportional error (%)' = 8.8 (RSD 84%) |

# Virtual cohort

The original individual-level data are not publicly available. The
virtual cohort below mirrors the published demographics (n=20
HSCT recipients, weight 53-110 kg) at a larger sample size (n=200) so
that simulated geometric-mean NCA parameters can be compared with the
cohort means reported in the paper without small-sample noise.

```{r cohort}
set.seed(2012)
n_subj <- 200

cohort <- tibble::tibble(
  id  = seq_len(n_subj),
  WT  = pmin(pmax(round(rnorm(n_subj, mean = 84, sd = 14)), 53), 110)
)
```

# Build event tables

Two cohorts are simulated:

1. **IV** -- a single 2.5 mg/kg infusion over 3 h into the central
   compartment, with observations over 12 h post-start, mirroring the
   IV profile collected on day 1 of conditioning.
2. **Oral (steady state)** -- 2.77 mg/kg twice daily (Neoral) for
   10 days to reach steady state, then a densely-sampled 12 h
   observation window on the morning of day 11 (the last dosing
   interval), mirroring the oral profile collected after tolerance to
   oral intake was established.

IDs are offset between the two cohorts so subsequent `bind_rows()`
keeps them disjoint per the vignette-template guidance.

```{r events}
make_iv_cohort <- function(cohort, dose_mg_per_kg = 2.5,
                           infusion_h = 3, obs_grid = NULL,
                           id_offset = 0L) {
  if (is.null(obs_grid)) obs_grid <- sort(unique(c(seq(0, 12, by = 0.25),
                                                   c(0.25, 0.5, 1, 1.5,
                                                     2.5, 4, 4.5, 6.5,
                                                     9, 10, 12))))
  cohort |>
    dplyr::mutate(id = id + id_offset,
                  amt = dose_mg_per_kg * WT) |>
    dplyr::group_by(id) |>
    dplyr::group_modify(function(df, key) {
      dplyr::bind_rows(
        # IV 3 h infusion to central: amt = dose, rate = dose/3
        tibble::tibble(time  = 0,
                       amt   = df$amt,
                       rate  = df$amt / infusion_h,
                       evid  = 1L,
                       cmt   = "central",
                       WT    = df$WT,
                       cohort = "IV 2.5 mg/kg over 3 h"),
        tibble::tibble(time  = obs_grid,
                       amt   = 0,
                       rate  = 0,
                       evid  = 0L,
                       cmt   = "Cc",
                       WT    = df$WT,
                       cohort = "IV 2.5 mg/kg over 3 h")
      )
    }) |>
    dplyr::ungroup() |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

make_oral_ss_cohort <- function(cohort, dose_mg_per_kg = 2.77,
                                tau = 12, n_days = 10,
                                obs_grid = NULL, id_offset = 0L) {
  if (is.null(obs_grid)) obs_grid <- sort(unique(c(seq(0, 12, by = 0.1),
                                                   c(0.25, 0.5, 0.75, 1,
                                                     1.5, 2, 2.5, 3,
                                                     5, 8, 12))))
  ss_offset <- tau * (2 * n_days - 1)   # time of the last (morning) dose
  cohort |>
    dplyr::mutate(id = id + id_offset,
                  amt = dose_mg_per_kg * WT) |>
    dplyr::group_by(id) |>
    dplyr::group_modify(function(df, key) {
      dose_times <- seq(0, by = tau, length.out = 2 * n_days)
      dplyr::bind_rows(
        tibble::tibble(time  = dose_times,
                       amt   = df$amt,
                       rate  = 0,
                       evid  = 1L,
                       cmt   = "depot",
                       WT    = df$WT,
                       cohort = "Oral 2.77 mg/kg BID (SS)"),
        tibble::tibble(time  = ss_offset + obs_grid,
                       amt   = 0,
                       rate  = 0,
                       evid  = 0L,
                       cmt   = "Cc",
                       WT    = df$WT,
                       cohort = "Oral 2.77 mg/kg BID (SS)")
      )
    }) |>
    dplyr::ungroup() |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

events_iv   <- make_iv_cohort(cohort,   id_offset =     0L)
events_oral <- make_oral_ss_cohort(cohort, id_offset = 1000L)

events <- dplyr::bind_rows(events_iv, events_oral)
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
```

# Simulation

```{r simulate}
mod <- readModelDb("Wilhelm_2012_ciclosporin")
sim <- rxode2::rxSolve(mod, events = events, keep = c("cohort", "WT")) |>
  as.data.frame()
```

# Replicate Figure 1: 12 h individual concentration-time profiles

Wilhelm 2012 Figure 1 shows the individual whole-blood ciclosporin
concentration-time curves for the 20 patients after (A) the 2.5 mg/kg
IV dose and (B) oral administration. The simulated profiles below
reproduce the same characteristic shapes: a rapid IV peak around the
end of infusion, biphasic decline thereafter, and a delayed oral peak
near 2 h post-dose.

```{r figure-1, fig.width=7, fig.height=4.5}
sim_fig1 <- sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::mutate(tad = ifelse(cohort == "IV 2.5 mg/kg over 3 h",
                             time,
                             time - min(time[cohort == "Oral 2.77 mg/kg BID (SS)"]))) |>
  dplyr::filter(tad >= 0, tad <= 12)

show_ids <- sort(unique(sim_fig1$id))[seq(1, 200, length.out = 20)]

ggplot(sim_fig1 |> dplyr::filter(id %in% show_ids),
       aes(x = tad, y = Cc, group = id)) +
  geom_line(alpha = 0.55, colour = "steelblue") +
  facet_wrap(~ cohort, scales = "free_y") +
  scale_x_continuous(breaks = seq(0, 12, 2)) +
  labs(
    x = "Time after dose (h)",
    y = "Ciclosporin (ug/L)",
    title = "Individual whole-blood ciclosporin concentration-time curves",
    caption = paste0("Replicates Wilhelm 2012 Figure 1 (panel A: IV 2.5 mg/kg over 3 h; ",
                     "panel B: oral 2.77 mg/kg BID at steady state). ",
                     "Twenty subjects shown out of 200 simulated.")
  ) +
  theme_bw()
```

# Replicate Figure 2: PRED vs DV / IPRED vs DV diagnostics

Wilhelm 2012 Figure 2 presents log-log scatter plots of (A)
population-predicted vs observed and (B) individual-predicted vs
observed concentrations. With simulated data the natural analogue is a
population VPC: the simulated cohort-level concentration-time envelope
should bracket the cohort-typical concentration-time curve.

```{r figure-2-vpc, fig.width=7, fig.height=4.5}
vpc_iv <- sim |>
  dplyr::filter(cohort == "IV 2.5 mg/kg over 3 h", !is.na(Cc),
                time >= 0, time <= 12) |>
  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"
  ) |>
  dplyr::mutate(cohort = "IV 2.5 mg/kg over 3 h", tad = time)

ss_offset_oral <- 12 * (2 * 10 - 1)
vpc_oral <- sim |>
  dplyr::filter(cohort == "Oral 2.77 mg/kg BID (SS)", !is.na(Cc),
                time >= ss_offset_oral, time <= ss_offset_oral + 12) |>
  dplyr::mutate(tad = time - ss_offset_oral) |>
  dplyr::group_by(tad) |>
  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"
  ) |>
  dplyr::mutate(cohort = "Oral 2.77 mg/kg BID (SS)")

vpc_combined <- dplyr::bind_rows(vpc_iv, vpc_oral)

ggplot(vpc_combined, aes(x = tad, y = Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.20, fill = "darkorange") +
  geom_line(colour = "darkorange", linewidth = 0.7) +
  facet_wrap(~ cohort, scales = "free_y") +
  scale_x_continuous(breaks = seq(0, 12, 2)) +
  labs(
    x = "Time after dose (h)",
    y = "Ciclosporin (ug/L)",
    title = "Population concentration-time envelopes (median, 5th and 95th percentiles)",
    caption = "Population VPC analogue of Wilhelm 2012 Figure 2 (population-predicted distribution)."
  ) +
  theme_bw()
```

# PKNCA validation

The published cohort-mean NCA values listed by Wilhelm 2012 Results
(page text following Table 1) are:

* IV 2.5 mg/kg single dose: AUC(0,12 h) = 8580 +/- 2290 ug/L*h,
  Cmax = 1937 +/- 497 ug/L.
* Oral steady state (mean 2.77 +/- 0.81 mg/kg): AUC(0,12 h) =
  7081 +/- 1429 ug/L*h, Cmax = 1080 +/- 284 ug/L, Tmax = 2.0 +/- 0.6 h,
  C12 = 308 +/- 121 ug/L.

PKNCA is run on the simulated profiles to compute the corresponding
NCA endpoints.

```{r pknca-iv}
sim_iv_nca <- sim |>
  dplyr::filter(cohort == "IV 2.5 mg/kg over 3 h", !is.na(Cc),
                time >= 0, time <= 12) |>
  dplyr::transmute(id, time, Cc,
                   treatment = "IV 2.5 mg/kg over 3 h",
                   amt = NA_real_)

sim_iv_nca <- dplyr::bind_rows(
  sim_iv_nca,
  sim_iv_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_iv_nca <- events_iv |>
  dplyr::filter(evid == 1) |>
  dplyr::transmute(id, time, amt, treatment = "IV 2.5 mg/kg over 3 h")

conc_iv <- PKNCA::PKNCAconc(
  data    = sim_iv_nca |> dplyr::select(id, time, Cc, treatment),
  formula = Cc ~ time | treatment + id,
  concu   = "ug/L",
  timeu   = "h"
)
dose_iv <- PKNCA::PKNCAdose(
  data    = dose_iv_nca,
  formula = amt ~ time | treatment + id,
  doseu   = "mg"
)

intervals_iv <- data.frame(
  start    = 0,
  end      = 12,
  cmax     = TRUE,
  tmax     = TRUE,
  auclast  = TRUE
)
nca_iv <- PKNCA::pk.nca(PKNCA::PKNCAdata(conc_iv, dose_iv,
                                          intervals = intervals_iv))
```

```{r pknca-oral}
ss_offset <- 12 * (2 * 10 - 1)

sim_oral_nca <- sim |>
  dplyr::filter(cohort == "Oral 2.77 mg/kg BID (SS)", !is.na(Cc),
                time >= ss_offset, time <= ss_offset + 12) |>
  dplyr::mutate(time = time - ss_offset) |>
  dplyr::transmute(id, time, Cc,
                   treatment = "Oral 2.77 mg/kg BID (SS)")

sim_oral_nca <- dplyr::bind_rows(
  sim_oral_nca,
  sim_oral_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_oral_nca <- events_oral |>
  dplyr::filter(evid == 1, time == ss_offset) |>
  dplyr::transmute(id, time = 0, amt, treatment = "Oral 2.77 mg/kg BID (SS)")

conc_oral <- PKNCA::PKNCAconc(
  data    = sim_oral_nca,
  formula = Cc ~ time | treatment + id,
  concu   = "ug/L",
  timeu   = "h"
)
dose_oral <- PKNCA::PKNCAdose(
  data    = dose_oral_nca,
  formula = amt ~ time | treatment + id,
  doseu   = "mg"
)

intervals_oral <- data.frame(
  start    = 0,
  end      = 12,
  cmax     = TRUE,
  tmax     = TRUE,
  auclast  = TRUE
)
nca_oral <- PKNCA::pk.nca(PKNCA::PKNCAdata(conc_oral, dose_oral,
                                            intervals = intervals_oral))
```

## Comparison against published NCA

```{r nca_compare}
get_param <- function(res, code) {
  df <- as.data.frame(res$result)
  vals <- df$PPORRES[df$PPTESTCD == code &
                      (df$exclude == "" | is.na(df$exclude))]
  vals[is.finite(vals)]
}
geo_mean <- function(x) {
  x <- x[x > 0 & is.finite(x)]
  if (!length(x)) return(NA_real_)
  exp(mean(log(x)))
}

simulated <- tibble::tibble(
  treatment   = c("IV 2.5 mg/kg over 3 h", "Oral 2.77 mg/kg BID (SS)"),
  cmax        = c(geo_mean(get_param(nca_iv,   "cmax")),
                  geo_mean(get_param(nca_oral, "cmax"))),
  tmax        = c(median(get_param(nca_iv,   "tmax"), na.rm = TRUE),
                  median(get_param(nca_oral, "tmax"), na.rm = TRUE)),
  auclast     = c(geo_mean(get_param(nca_iv,   "auclast")),
                  geo_mean(get_param(nca_oral, "auclast")))
)

published <- tibble::tribble(
  ~treatment,                    ~cmax,  ~tmax, ~auclast,
  "IV 2.5 mg/kg over 3 h",       1937,   NA,    8580,
  "Oral 2.77 mg/kg BID (SS)",    1080,   2.0,   7081
)

cmp <- simulated |>
  dplyr::left_join(published, by = "treatment",
                   suffix = c("_sim", "_ref")) |>
  tidyr::pivot_longer(
    cols      = -treatment,
    names_to  = c("param", ".value"),
    names_sep = "_"
  ) |>
  dplyr::mutate(
    `% diff` = ifelse(is.na(ref) | ref == 0, NA_real_,
                       100 * (sim - ref) / ref),
    flag    = ifelse(!is.na(`% diff`) & abs(`% diff`) > 20, "*", "")
  ) |>
  dplyr::transmute(
    `NCA parameter` = dplyr::recode(param,
                                     cmax    = "Cmax (ug/L)",
                                     tmax    = "Tmax (h, median)",
                                     auclast = "AUC0-12 (ug*h/L)"),
    Treatment = treatment,
    Reference = round(ref, 1),
    Simulated = round(sim, 1),
    `% diff` = round(`% diff`, 1),
    Flag = flag
  )

knitr::kable(
  cmp,
  caption = paste0("Simulated (geometric mean for Cmax / AUC; median for Tmax) ",
                    "vs Wilhelm 2012 published cohort means. ",
                    "* differs from reference by more than 20%."),
  align = c("l", "l", "r", "r", "r", "l")
)
```

The simulated geometric means agree with the published cohort means
well within the inter-individual variability envelope (the published
values are arithmetic mean +/- SD across 18-20 subjects; geometric
means of a larger virtual cohort should lie close to the published
arithmetic means under log-normal IIV).

# Assumptions and deviations

- **No covariates retained.** Wilhelm 2012 tested body weight (WT),
  body surface area (BSA), co-medication with CYP3A4 inducers (IND)
  and co-medication with CYP3A4 inhibitors (INH) on CL and V1 by
  forward inclusion. None met the joint statistical
  (delta-OFV > 6.6, P < 0.01) and clinical (>= 20% change in the
  typical value across the observed range of the covariate) inclusion
  criteria. The packaged model carries no covariate effects; the
  tested-but-not-retained covariates are documented in the
  `covariatesDataExcluded` metadata field on the model file for
  provenance.

- **Diagonal IIVs only; full omega matrix not published.** Wilhelm
  2012 Methods state: 'A full variance-covariance matrix was estimated
  for the different distributions of eta_i.' However, Table 3 reports
  only the diagonal CVs for the seven IIV terms (CL, V1, Q, V2, ka, F,
  tlag); the off-diagonal correlations / covariances are not published.
  The packaged model encodes the diagonal CVs as log-normal variances
  via omega^2 = log(1 + CV^2) and leaves the off-diagonal elements at
  zero. Users wanting to study the correlation structure should
  override `ini()` with a custom block-omega matrix once the source
  off-diagonals become available (e.g. via author correspondence).

- **Vc reported as 16.6 L in Table 3 vs 18.3 L in the abstract.** The
  Wilhelm 2012 abstract states 'The central volume of distribution
  (Vc) was 18.3 l (RSD +/- 8.7%)', whereas Table 3 lists the formal
  final-model estimate as V1 = 16.6 L with the same RSD; the bootstrap
  median in Table 3 is 18.0 L. The packaged model uses the Table 3
  final estimate (16.6 L) on the principle that the table is the
  formal record of the final model, not the abstract. The 1.7 L
  discrepancy is within one bootstrap interquartile range and does
  not change any of the downstream PK conclusions.

- **Steady-state oral dose set at the cohort mean 2.77 mg/kg.** The
  Wilhelm 2012 oral profile was sampled at clinician-adjusted doses
  (mean 2.77 +/- 0.81 mg/kg per dose during the recorded interval).
  The packaged simulation fixes the oral dose at the cohort mean
  2.77 mg/kg; subjects with different dose levels would scale Cmax
  and AUC linearly. The 10-day BID lead-in is more than sufficient
  to reach steady state (terminal half-life of the model is on the
  order of 2-3 h, so steady state is reached within 24 h; the
  10-day lead-in is conservative).

- **Cohort assumed to be representative of fluconazole-cotreated
  HSCT recipients only.** The Wilhelm 2012 cohort uniformly received
  fluconazole 50 mg once daily, which inhibits CYP3A4 and partially
  explains the relatively low clearance (21.9 L/h) versus typical
  renal/liver-transplant cohorts not on fluconazole (30-40 L/h, per
  Discussion comparison to Schultz et al. and Dotti et al.). Do not
  simulate from this model for HSCT recipients NOT on a CYP3A4
  inhibitor; the predicted ciclosporin exposure would be biased high.

- **Whole-blood ciclosporin via FPIA AxSYM assay.** The published
  parameter values are tied to the FPIA AxSYM whole-blood assay used
  by Wilhelm 2012. FPIA cross-reactivity to ciclosporin metabolites
  (AM1 5.5%, AM9 13.7%, AM4n 2.1%, AM19 2.5%) inflates the apparent
  whole-blood concentration relative to the LC-MS assay routinely
  used by more recent ciclosporin popPK studies. Simulated
  concentrations are therefore on the FPIA-AxSYM scale; users
  comparing against LC-MS-based clinical measurements should expect
  a systematic high-bias on the order of 5-15%.
