---
title: "Savic_2017_cladribine"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Savic_2017_cladribine}
  %\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("Savic_2017_cladribine"))$reference`
* Description: `r rxode2::rxode(readModelDb("Savic_2017_cladribine"))$description`
* Article: <https://doi.org/10.1007/s40262-017-0516-6>

# Population

The Savic 2017 dataset pools four clinical studies in patients with
relapsing-remitting multiple sclerosis (Table 2 of the source paper):

* Study 25803 (phase I, n = 16): randomized two-period crossover; 3 mg IV
  infusion over 1 h or 10 mg single oral tablet, with parent cladribine (CdA)
  and metabolite 2-chloroadenine (CAde) plasma and urine sampling.
* Study 26127 (phase I, n = 16): food-effect crossover, single 10 mg oral
  tablet in fed vs fasted sequence.
* Study 26486 (phase I, n = 16): cladribine-and-interferon-beta-1a drug-drug
  interaction study; oral cladribine 10 mg daily for 5 consecutive days
  with subcutaneous IFN beta-1a 8.8/22/44 ug three times weekly.
* Study 25643 (CLARITY phase III, n = 125; NCT00213135): 96-week,
  double-blind, placebo-controlled trial of cumulative 3.5 or 5.25 mg/kg
  oral cladribine tablets over 2 years.

Pooled demographics (source paper Table 3): N = 173, age range
19-65 years (median 40), body weight 48.5-116.1 kg (median 69.2),
65.9 % female, Cockcroft-Gault creatinine clearance 49.6-244.4 mL/min
(median 107.9). The original dataset contained 4790 concentration-time
records; 45 % were excluded (96.5 % of those for BLQ), leaving 2619
plasma and urine records for the population analysis.

```{r population-metadata}
str(rxode2::rxode(readModelDb("Savic_2017_cladribine"))$population)
```

# Source trace

Per-parameter origin is recorded inline next to each `ini()` entry in
`inst/modeldb/specificDrugs/Savic_2017_cladribine.R`. Consolidated here for
review:

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lka` (Ka, fasted)                                  | 1.08 h^-1                       | Table 4 (RSE 21.14 %) |
| `e_fed_ka` (FED on log-Ka)                          | log(1.03 / 1.08) = -0.0474      | Table 4 ("unknown/fed" Ka = 1.03 h^-1) |
| `logitf` (logit of F, fasted)                       | logit(0.456) = -0.1765          | Table 4 (RSE 7.03 %); footnote d "Variance on a logit scale" |
| `e_fed_logitf` (FED on logit-F)                     | logit(0.4) - logit(0.456) = -0.229 | Table 4 ("unknown/fed" F = 0.4) |
| `lmtt_fed` (lag time approximation, fed)            | log(0.910 h) = -0.0943          | Table 4 mean transit time (fed state) = 0.910 h |
| `lcl_renal_typ` (CLR at typical CRCL)               | log(22.2) = 3.100               | Table 4 CLR coefficient 3.52 dimensionless; CLR = 22.2 L/h at CLCR = 6.31 L/h = 105.2 mL/min |
| `lcl_nonrenal` (CLNR, monotherapy)                  | log(23.4) = 3.153               | Table 4 (RSE 9.58 %) |
| `e_ifn_clnr` (IFN beta-1a fractional CLNR increase) | 0.21                            | Table 4 ("Fold increase in CLNR in the presence of IFN beta-1a" = 1.21) |
| `lvc` (V1)                                          | log(44.0) = 3.784               | Table 4 (RSE 22.77 %) |
| `lq` (Q3)                                           | log(14.3) = 2.660               | Table 4 (RSE 7.73 %) |
| `lq2` (Q4)                                          | log(53.7) = 3.983               | Table 4 (RSE 19.06 %) |
| `lvp` (V3)                                          | log(347) = 5.849                | Table 4 (RSE 6.07 %) |
| `lvp2` (V4)                                         | log(89.5) = 4.494               | Table 4 (RSE 7.97 %) |
| `etalka` (BSV Ka, variance)                         | 0.102                           | Table 4 (RSE 62.53 %) |
| `etalogitf` (BSV F, logit-scale variance)           | 0.223                           | Table 4 (RSE 19.13 %, footnote d) |
| `etalcl_nonrenal` (BSV CLNR, variance)              | 0.00574                         | Table 4 (RSE 64.41 %) |
| `etalvc` (BSV V1, variance)                         | 0.209                           | Table 4 (RSE 72.51 %) |
| `etalq` (BSV Q3/Q4/V3/V4 shared, variance)          | 0.0365                          | Table 4 (RSE 28.47 %, single shared eta) |
| `propSd` (proportional residual, plasma)            | 0.353                           | Table 4 RUV plasma (oral, study 25643) = 35.3 % |
| Three-compartment ODE structure                     | n/a                             | Methods 2.3.1 ("Model I A three-compartment model best described the PK of CdA") and Figure 1 |
| CLR = coefficient * CRCL                            | n/a                             | Methods 3.2.1 and Table 4 footnotes b, c |
| IFN beta-1a effect on CLNR                          | multiplicative *1.21            | Methods 3.2.1 (Model IV) |

# Virtual cohort

Original observed data are not publicly available. The figures below use
a virtual population whose covariate distributions approximate the
published pooled demographics (Table 3).

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

n_sub <- 100L

make_cohort <- function(n, label, fed, ifnb1a, crcl_mean = 107.9,
                        crcl_sd = 35, id_offset = 0L) {
  ids <- id_offset + seq_len(n)
  base <- tibble::tibble(
    id            = ids,
    cohort        = label,
    CRCL          = pmax(40, rnorm(n, mean = crcl_mean, sd = crcl_sd)),
    FED           = as.integer(fed),
    CONMED_IFNB1A = as.integer(ifnb1a)
  )
  obs_grid <- c(0, 0.25, 0.5, 1, 1.5, 2, 3, 4, 6, 8, 12, 16, 24, 36, 48, 72)
  expand_one <- function(row) {
    dose <- data.frame(id = row$id, cohort = row$cohort, time = 0,
                       evid = 1L, amt = 10, cmt = "depot",
                       CRCL = row$CRCL, FED = row$FED,
                       CONMED_IFNB1A = row$CONMED_IFNB1A)
    obs  <- data.frame(id = row$id, cohort = row$cohort, time = obs_grid,
                       evid = 0L, amt = 0, cmt = "Cc",
                       CRCL = row$CRCL, FED = row$FED,
                       CONMED_IFNB1A = row$CONMED_IFNB1A)
    rbind(dose, obs)
  }
  rows <- split(base, base$id)
  do.call(rbind, lapply(rows, expand_one))
}

events <- dplyr::bind_rows(
  make_cohort(n_sub, "10 mg oral, fasted, no IFN",       fed = 0, ifnb1a = 0,
              id_offset =      0L),
  make_cohort(n_sub, "10 mg oral, fed, no IFN",          fed = 1, ifnb1a = 0,
              id_offset = 1*n_sub),
  make_cohort(n_sub, "10 mg oral, fasted, + IFN beta-1a", fed = 0, ifnb1a = 1,
              id_offset = 2*n_sub)
)
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
```

# Simulation

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

# Replicate published figures

## Figure 2 / Figure 4 -- typical-value plasma concentration profile

Savic 2017 Figure 2c shows the VPC for a single oral 10 mg dose in fasted
state. We reproduce the typical-value (median) cladribine plasma
concentration over 72 hours.

```{r figure-2c-typical}
typical_summary <- sim |>
  dplyr::filter(time > 0) |>
  dplyr::group_by(cohort, time) |>
  dplyr::summarise(
    median_Cc = median(Cc),
    p05_Cc    = quantile(Cc, 0.05),
    p95_Cc    = quantile(Cc, 0.95),
    .groups   = "drop"
  )

ggplot(typical_summary, aes(time, median_Cc, colour = cohort, fill = cohort)) +
  geom_ribbon(aes(ymin = p05_Cc, ymax = p95_Cc), alpha = 0.2, colour = NA) +
  geom_line(linewidth = 0.8) +
  scale_y_log10() +
  scale_x_continuous(breaks = c(0, 6, 12, 24, 36, 48, 72)) +
  labs(x = "Time after dose (h)", y = "Cladribine plasma conc (ng/mL)",
       title = "Cladribine 10 mg oral, simulated VPC by cohort",
       subtitle = "Replicates the shape of Savic 2017 Figure 2 / Figure 4",
       colour = NULL, fill = NULL) +
  theme_minimal() +
  theme(legend.position = "bottom")
```

## Renal-impairment exposure scaling

The Discussion of Savic 2017 reports model-predicted decreases in total
cladribine clearance of 19 %, 30 %, and 40 % for CLCR = 65, 40, and
20 mL/min, respectively. Reproduce by simulating a 10 mg oral fasted
dose across a CLCR grid.

```{r renal-impairment}
ri_grid <- c(20, 40, 65, 105.2)
ri_events <- dplyr::bind_rows(lapply(seq_along(ri_grid), function(i) {
  crcl <- ri_grid[i]
  id_off <- (i - 1) * 1L
  obs_grid <- c(seq(0, 24, by = 0.25), seq(25, 168, by = 1))
  data.frame(id    = c(id_off + 1L, rep(id_off + 1L, length(obs_grid))),
             time  = c(0, obs_grid),
             evid  = c(1L, rep(0L, length(obs_grid))),
             amt   = c(10, rep(0, length(obs_grid))),
             cmt   = c("depot", rep("Cc", length(obs_grid))),
             CRCL  = crcl,
             FED   = 0L,
             CONMED_IFNB1A = 0L,
             crcl_label = paste0("CLCR = ", crcl, " mL/min"))
}))

mod_typ <- rxode2::zeroRe(mod())
ri_sim <- as.data.frame(rxode2::rxSolve(mod_typ, events = ri_events,
                                         keep = c("crcl_label")))

ri_total_cl <- ri_sim |>
  dplyr::filter(time == 0.25) |>           # any post-dose time works at typical values
  dplyr::transmute(crcl_label, total_cl = cl) |>
  dplyr::distinct()

baseline_cl <- ri_total_cl$total_cl[ri_total_cl$crcl_label == "CLCR = 105.2 mL/min"]
ri_total_cl$pct_decrease <- round(
  100 * (baseline_cl - ri_total_cl$total_cl) / baseline_cl, 1)
knitr::kable(ri_total_cl,
             caption = "Predicted total CL vs published Discussion targets (Savic 2017): 19 % at 65, 30 % at 40, 40 % at 20 mL/min.")

ggplot(ri_sim, aes(time, Cc, colour = crcl_label)) +
  geom_line(linewidth = 0.8) +
  scale_y_log10() +
  scale_x_continuous(breaks = c(0, 6, 12, 24, 48, 72, 96, 120, 168)) +
  labs(x = "Time after dose (h)", y = "Cladribine plasma conc (ng/mL)",
       title = "Renal-impairment exposure scaling (typical values, 10 mg oral, fasted)",
       colour = NULL) +
  theme_minimal() +
  theme(legend.position = "bottom")
```

# PKNCA validation

Use PKNCA for Cmax, Tmax, AUClast, and AUCinf on the three simulated
cohorts (fasted monotherapy, fed monotherapy, fasted + IFN beta-1a). Per
the source paper Table 4 footnote a, the typical-patient F and Ka used
in the model imply published AUC and Cmax values consistent with the
sparse-NCA values reported in the literature for cladribine 10 mg oral.

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

dose_df <- events |>
  dplyr::filter(evid == 1) |>
  dplyr::select(id, time, amt, cohort)

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

intervals <- data.frame(
  start      = 0,
  end        = Inf,
  cmax       = TRUE,
  tmax       = TRUE,
  auclast    = 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_summary <- summary(nca_res)
knitr::kable(nca_summary,
             caption = "Simulated NCA parameters by cohort (Savic 2017 cladribine 10 mg oral).")
```

## Comparison against published model-implied targets

Savic 2017 does not publish a per-cohort NCA table; the paper validates
the model via VPCs in Figures 2-4. Three internally-consistent checks
are available:

* Bioavailability ratio fed vs fasted: F(fed) / F(fasted) = 0.4 / 0.456
  ~= 0.877 (paper says food reduced bioavailability by ~12 %; the model
  reproduces this ratio).
* IFN beta-1a effect: total CL increases by ~11 % when CONMED_IFNB1A = 1
  ((22.2 + 23.4 * 1.21) / (22.2 + 23.4) ~= 1.108; paper Abstract: "an
  increase of 11.1 % in total clearance").
* Renal-impairment scaling: total CL decreases by ~30 % at CLCR =
  40 mL/min (paper Discussion: 30 %). See the renal-impairment table
  above.

# Assumptions and deviations

Several simplifications were made when packaging the Savic 2017 final
model. The full source-paper structure can be recovered by extending
the packaged model where needed.

* **Parent-only model.** Savic 2017 jointly fits cladribine and its
  primary metabolite 2-chloroadenine (CAde). The packaged model
  implements the cladribine parent kinetics only. The CAde compartment
  was omitted because Table 4 reports CAde apparent CL/F_m and V/F_m
  *corrected for the fraction metabolized* (footnote f), and F_m is
  not separately identifiable in the source data. Without F_m the CAde
  compartment cannot be simulated on an absolute concentration scale.
  Users needing metabolite predictions should extend this model by
  adding a CAde compartment with formation rate F_m * CLNR(CdA) /
  Vc(CdA) and the apparent Vmax/Km saturable elimination from Table 4,
  acknowledging the F_m ambiguity.
* **Transit-compartment fed-state absorption simplified to a lag time.**
  Savic 2017 best fit the fed-state absorption with a transit-compartment
  chain (mean transit time MTT = 0.910 h, n = 2.24 compartments). The
  packaged model approximates this with a fed-state-only lag time
  alag(depot) = MTT * FED. The approximation preserves the absorption-
  delay magnitude but loses the gamma-distributed input-rate shape; the
  source paper notes the transit model improved fit over a pure lag-time
  approach, so for downstream simulation needing the precise absorption
  profile users should override with `transit(2.24, 0.910)` gated on
  FED.
* **Phase III-specific 0.319 h lag time omitted.** Table 4 carries a
  separate 0.319 h lag time fitted only on the CLARITY (study 25643)
  data. This is a study-specific artefact (different sampling /
  dosing logistics in the larger phase III) rather than a
  generalizable patient-level effect; it is not reproduced here.
* **Single proportional residual error magnitude.** Table 4 reports
  four residual-error magnitudes by route/study (IV 20.0 %, oral
  studies 25803 34.7 %, 26127 + 26486 22.8 %, 25643 35.3 %) plus a
  separate BSV of 0.159 on the residual error. The packaged model
  uses the phase III value (35.3 %) as the largest cohort and most
  conservative simulation choice, and omits the BSV-on-residual term.
* **Single shared eta on distribution parameters.** Table 4 reports a
  single BSV of 0.0365 for the four parameters Q3, Q4, V3, V4, which is
  consistent with a single NONMEM ETA applied multiplicatively to all
  four. The packaged model implements this via one `etalq` random
  effect added to each of `lq`, `lq2`, `lvp`, and `lvp2`. Users who
  prefer four diagonal etas (e.g. for sensitivity analysis) can swap
  the shared eta for `etalq ~ 0.0365`, `etalq2 ~ 0.0365`,
  `etalvp ~ 0.0365`, `etalvp2 ~ 0.0365` and re-load the model.
* **Logit-scale BSV on bioavailability.** Savic 2017 Table 4 footnote
  d notes "Variance on a logit scale" for the BSV of F. The packaged
  model implements this by adding `etalogitf` inside the logit
  transform, then applying the inverse logit to keep F bounded in
  (0, 1). The variance 0.223 is on the logit scale, not the linear
  F scale.
* **"Unknown" food-state records are pooled with fed records.** Table 4
  reports a single "unknown/fed state" column for Ka and F. Records in
  studies 26486 and 25643 that did not record food state share
  parameters with the fed state. In the packaged model, FED = 1 covers
  both confirmed-fed and unknown-food-state records (see the
  `covariateData[[FED]]$notes` field).
* **CRCL canonical column.** The source paper's CLCR (raw
  Cockcroft-Gault, mL/min, not BSA-normalized) is stored under the
  canonical CRCL column per `inst/references/covariate-columns.md`
  (CRCL accepts raw Cockcroft-Gault mL/min when the paper does not
  apply BSA normalization). The model uses CRCL directly with
  reference value 105.2 mL/min (the paper's typical patient CLCR =
  6.31 L/h).
* **CONMED_IFNB1A canonical covariate.** Adding this model also
  introduces the new canonical column CONMED_IFNB1A
  (concomitant subcutaneous IFN beta-1a coadministration; binary,
  general scope, reference 0 = no concomitant IFN). The register
  entry is committed alongside this model in
  `inst/references/covariate-columns.md`.
