---
title: "Cyclosporine (Philippe 2015)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Cyclosporine (Philippe 2015)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4)
library(nlmixr2lib)
library(PKNCA)
library(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
```

# Model and source

* Citation: `r readModelDb("Philippe_2015_cyclosporine")()$reference`
* Description: `r readModelDb("Philippe_2015_cyclosporine")()$description`
* Article: <https://doi.org/10.1208/s12248-015-9779-8>

The packaged model is a pediatric PK-PD-time-to-event model for oral
cyclosporine (CsA) in children with severe aplastic anemia (SAA). The PK
is a two-compartment model with first-order absorption and an absorption
lag time, with absorption parameters (F, Tlag, ka) fixed from the
literature and body weight scaling on V1, V2, Cl, and Q. The
pharmacodynamics combines an interface model (Eq. 5 of the paper) that
filters the trough concentration through a window-of-efficacy production
term with a time-to-event hazard (Eq. 6) for the neutrophil response.

# Population

The source population was 21 pediatric patients (23 treatment courses,
aged 1-15 years, weight 9.8-79.3 kg) with severe aplastic anemia treated
in Lyon (n = 22 courses) and Grenoble (n = 1 course), France, between
1998 and 2013. Three patients who experienced first-course failure were
modelled as independent individuals at the second course. The
immunosuppressive regimen was anti-thymocyte globulin (160 mg/kg horse
ATG or 18.75 mg/kg rabbit ATG) plus oral cyclosporine started at 5
mg/kg every 12 hours, subsequently adjusted to an initial trough target
of 150 ng/mL. 341 trough concentrations were used for the PK model
(median 14 per patient). 15/23 courses (65.2%) achieved a neutrophil
response (ANC > 0.5e9/L on two consecutive measurements) with a mean
time-to-response of 69 days (range 19-182). Baseline demographics are
in Table I of the source paper; the same fields are available
programmatically via `readModelDb("Philippe_2015_cyclosporine")$population`.

# Source trace

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

| Equation / parameter             | Value                | Source location                          |
|----------------------------------|----------------------|------------------------------------------|
| F (fixed)                        | 0.386                | Table II                                 |
| Tlag (fixed)                     | 0.648 h              | Table II                                 |
| ka (fixed)                       | 0.829 1/h            | Table II                                 |
| Cl (at WT = 34 kg)               | 7.2 L/h              | Table II                                 |
| V1 (at WT = 34 kg)               | 37.5 L               | Table II                                 |
| Q (at WT = 34 kg)                | 2.7 L/h              | Table II                                 |
| V2 (at WT = 34 kg)               | 1690 L               | Table II                                 |
| Allometric exponent on Cl, Q     | 0.75 (fixed)         | Methods Eq. 4                            |
| Allometric exponent on V1, V2    | 1 (fixed)            | Methods Eq. 4                            |
| Reference weight (BWpop)         | 34 kg                | Table I (cohort mean)                    |
| IIV Cl (CV%)                     | 31.4                 | Table II                                 |
| IIV V1 (CV%)                     | 30                   | Table II                                 |
| IIV Q  (CV%)                     | 66.2                 | Table II                                 |
| IIV V2 (CV%)                     | 30                   | Table II                                 |
| Additive residual SD             | 0.03 mg/L            | Table II                                 |
| Proportional residual SD         | 25%                  | Table II                                 |
| Interface ODE for Ce             | n/a                  | Methods Eq. 5                            |
| Lower bound gamma1               | 87 ng/mL             | Table II                                 |
| Upper bound gamma2               | 120 ng/mL            | Table II                                 |
| Interface rate alpha             | 0.028 1/day          | Table II (= 25-day half-life)            |
| Hazard equation lambda(t)        | n/a                  | Methods Eq. 6                            |
| Baseline hazard lambda0          | 0.0027 1/day         | Table II                                 |
| Hazard slope on Ce               | 17.2 mL/ng           | Table II                                 |
| Probability of response F(t)     | 1 - exp(-cumhaz)     | Methods Eq. 7                            |

# Virtual cohort

The original observed data are not publicly available. The figures
below use a small virtual cohort whose weight distribution approximates
the published trial demographics (Table I, n = 23 courses, mean weight
34 kg, range 9.8-79.3 kg).

```{r cohort}
set.seed(2015)
n_subj <- 30L
cohort <- tibble(
  id = seq_len(n_subj),
  # Spread weights across the published range (9.8-79.3 kg). Use a
  # log-uniform draw so the distribution is roughly cohort-like across
  # young and older children rather than weighted toward the larger ages.
  WT = exp(runif(n_subj, log(9.8), log(79.3)))
)
summary(cohort$WT)
```

# Single-individual PK profile

We first demonstrate the PK behavior in a single typical individual at
the published initial regimen (5 mg/kg every 12 hours) for the first 7
days, then a longer run at a maintenance dose chosen to bring the
predicted trough into the published efficacy window (87-120 ng/mL).

```{r pk-typical, fig.height = 4}
mod <- readModelDb("Philippe_2015_cyclosporine")

# Typical individual: WT = 30 kg, 5 mg/kg q12h x 7 days. Zero out the
# random effects to show the typical-value profile.
mod_typical <- rxode2::zeroRe(mod)

dose_typical <- 5 * 30
ev_typical <- rxode2::et(
  amt = dose_typical, ii = 12, until = 24 * 7, cmt = "depot"
)
ev_typical <- rxode2::et(ev_typical, seq(0, 24 * 7, by = 0.25))
ev_typical$WT <- 30

sim_typical <- rxode2::rxSolve(mod_typical, ev_typical)

ggplot(sim_typical, aes(time, 1000 * Cc)) +
  geom_line() +
  geom_hline(yintercept = c(87, 120), linetype = "dashed", colour = "red") +
  labs(
    x = "Time (h)",
    y = "Predicted whole-blood CsA concentration (ng/mL)",
    title = "Typical PK profile at 5 mg/kg q12h (WT = 30 kg)",
    caption = "Dashed red lines mark the published efficacy window gamma1 = 87 ng/mL and gamma2 = 120 ng/mL."
  )
```

At the initial 5 mg/kg q12h regimen, the typical-value trough exceeds
the published efficacy window of 87-120 ng/mL, which is consistent with
the paper's argument that the historical target of 200-400 ng/mL is too
high.

# Replicate Figure 4: probability of response over time

The paper's Figure 4 is a VPC of the probability of response (Eq. 7
in the source paper, F(t) = 1 - exp(-integrated hazard)) over time. We
simulate the cohort at a maintenance dose that brings most subjects'
trough into the 87-120 ng/mL window, then report the median and 5/95
percentile envelope of F(t) = 1 - sur across the cohort.

```{r figure-4, fig.height = 4}
# 1.0 mg/kg q12h for 6 months. Use 6-hourly sampling for the long
# horizon so the simulation finishes well under the 5-minute vignette
# gate while still resolving the survival curve and the interface
# dynamics (alpha = 0.028 1/day -> 25-day half-life).
ev_pop <- cohort |>
  rowwise() |>
  do({
    sub <- .
    ev <- rxode2::et(amt = 1.0 * sub$WT, ii = 12, until = 24 * 180, cmt = "depot")
    ev <- rxode2::et(ev, seq(0, 24 * 180, by = 6))
    ev <- as.data.frame(ev)
    ev$id <- sub$id
    ev$WT <- sub$WT
    ev
  }) |>
  ungroup() |>
  arrange(id, time, evid)

sim_pop <- rxode2::rxSolve(mod, ev_pop, keep = "WT") |>
  as.data.frame()
sim_pop$prob_response <- 1 - sim_pop$sur
sim_pop$day <- sim_pop$time / 24

vpc <- sim_pop |>
  group_by(day) |>
  summarise(
    Q05 = quantile(prob_response, 0.05),
    Q50 = quantile(prob_response, 0.50),
    Q95 = quantile(prob_response, 0.95),
    .groups = "drop"
  )

ggplot(vpc, aes(day, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25) +
  geom_line() +
  labs(
    x = "Time since treatment initiation (days)",
    y = "Probability of neutrophil response F(t) = 1 - sur",
    title = "Figure 4 analogue: probability of response over time",
    caption = "Replicates Figure 4 of Philippe 2015 (1.0 mg/kg q12h cohort, n = 30 simulated subjects)."
  )
```

The simulated median probability of response rises with the
characteristic delay of the interface model (effective-concentration
half-life of 25 days) and approaches ~60-70% by day 180, in line with
the paper's reported 65.2% response rate (Table I) and the median /
percentile pattern in Figure 4.

# PKNCA validation

The source paper does not report standard NCA parameters (Cmax, Tmax,
AUC, terminal half-life) but mentions that the cyclosporine blood
half-life is approximately 2.6 hours. We compute steady-state NCA from
the simulated cohort under the maintenance regimen above to confirm
the model produces PK behavior in that range.

```{r pknca}
# Use the last 12-hour dosing interval of the long simulation
# (days 6.5-7.0; covers t = 156 h to t = 168 h) so the PK has reached
# steady state. Resample at 0.5 h grid to give PKNCA enough resolution
# for Cmax / Tmax.
ev_nca <- cohort |>
  rowwise() |>
  do({
    sub <- .
    ev <- rxode2::et(amt = 1.0 * sub$WT, ii = 12, until = 24 * 7, cmt = "depot")
    ev <- rxode2::et(ev, seq(0, 24 * 7, by = 0.5))
    ev <- as.data.frame(ev)
    ev$id <- sub$id
    ev$WT <- sub$WT
    ev
  }) |>
  ungroup() |>
  arrange(id, time, evid)

sim_nca <- rxode2::rxSolve(mod, ev_nca, keep = "WT") |>
  as.data.frame() |>
  filter(!is.na(Cc), time >= 24 * 6.5, time <= 24 * 7)
sim_nca$conc_ngmL <- 1000 * sim_nca$Cc
sim_nca$treatment <- "1.0 mg/kg q12h"

conc_obj <- PKNCA::PKNCAconc(
  sim_nca, conc_ngmL ~ time | treatment + id
)

# Reconstruct the dose-time grid manually because rxode2::et collapses
# repeating doses into a single row with an ADDL count rather than one
# row per dose. PKNCA needs an explicit dose record for the start of
# the steady-state interval.
dose_df <- cohort |>
  mutate(
    time = 24 * 6,
    amt  = 1.0 * WT,
    treatment = "1.0 mg/kg q12h"
  ) |>
  select(id, time, amt, treatment)

dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id)

intervals <- data.frame(
  start            = 24 * 6,
  end              = 24 * 7,
  cmax             = TRUE,
  tmax             = TRUE,
  auclast          = TRUE,
  half.life        = TRUE
)

nca_data <- PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals)
nca_res  <- PKNCA::pk.nca(nca_data)
nca_tbl  <- as.data.frame(nca_res$result)

nca_tbl |>
  group_by(PPTESTCD) |>
  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"
  ) |>
  knitr::kable(
    caption = "Steady-state NCA summary across the simulated cohort (concentration in ng/mL, time in hours).",
    digits = 2
  )
```

The simulated terminal half-life is in the few-hours range that the
paper cites for CsA in blood. The dose was chosen so the typical
predicted trough sits near the lower end of the published efficacy
window (87 ng/mL).

# Assumptions and deviations

- **Ctrough vs continuous Cc.** The source paper's Eq. 5 uses the
  predicted trough concentration Ctrough (the model-predicted Cc value
  at the q12h trough times) as the input to the interface model. The
  packaged model uses the continuous predicted Cc (multiplied by 1000
  to convert mg/L to ng/mL) as a proxy for Ctrough, so a single ODE
  system can be simulated for arbitrary dosing schedules. This
  simplification is faithful in the limit of slow PD dynamics: the
  interface model's elimination half-life is approximately 25 days
  (alpha = 0.028 1/day), so PK oscillations on the 12-hour scale are
  heavily smoothed by the interface model. Users who require a strict
  trough-driven evaluation can simulate Cc continuously, sample at the
  trough times, then drive the interface ODE manually outside the
  packaged model.
- **Time axis in hours.** The source paper expresses the interface
  model and hazard equations in 1/day (alpha, lambda0). The packaged
  model uses time in hours to match the PK absorption parameters (ka,
  Tlag) so a single TIME axis works for both PK and PD. Inside the
  model, alpha and lambda0 are divided by 24 to give per-hour rates,
  and the sqrt production term in Eq. 5 is divided by 24 to convert
  from a per-day to a per-hour production rate.
- **Fixed absorption parameters.** F = 0.386, Tlag = 0.648 h, and
  ka = 0.829 1/h were fixed in the source paper from a prior PK study
  of cyclosporine in children (Methods page 1158, reference 24); the
  packaged model preserves this by wrapping each in `fixed()` in
  `ini()`.
- **Fixed allometric exponents.** The allometric exponents on Cl, Q
  (0.75) and on V1, V2 (1) are fixed at the standard physiological
  values per Methods Eq. 4 and wrapped in `fixed()` in `ini()`. The
  reference weight is the cohort mean of 34 kg (Table I).
- **No IIV on PD parameters.** The source paper reports "No IIV could
  be estimated on these parameters (IIV fixed to 0 for all the
  parameters)" for lambda0, alpha, slope, gamma1, gamma2 (Results
  page 1162). The packaged model omits any `eta` term on these
  parameters; their typical values are estimated (not fixed) and so
  are not wrapped in `fixed()`.
- **Cohort weight sampling for the vignette.** Table I reports the
  cohort weight range (9.8-79.3 kg) and mean (34 kg) but not the
  underlying distribution. The vignette draws weights from a
  log-uniform distribution over the reported range. This is a
  visualization convenience for the figure-4 analog; the model itself
  is independent of this choice.
- **Maintenance dose for the figure-4 analog.** The paper's protocol
  started at 5 mg/kg q12h and titrated to a 150 ng/mL trough target.
  Because individual titration histories are not published, the
  figure-4 analog instead uses a single fixed maintenance dose
  (1.0 mg/kg q12h) chosen so the typical-value trough sits near the
  lower end of the published efficacy window. The qualitative shape
  of the response-probability curve is what is being compared.
