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

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>"
)
```

```{r setup}
library(nlmixr2lib)
library(PKNCA)
library(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
```

# Model and source

```{r model-header, echo = FALSE, results = "asis"}
mod <- rxode2::rxode(readModelDb("Rosario_2015_vedolizumab"))
cat("* **Citation:** ", mod$reference, "\n", sep = "")
cat("* **Description:** ", mod$description, "\n", sep = "")
```

* **Article:** <https://doi.org/10.1111/apt.13243>
* **Corrigendum:** <https://doi.org/10.1111/apt.15571> — a text-only
  unit-typo fix (ng/mL → µg/mL) with no parameter-value changes.

# Population

The model was developed on 2700 evaluable PK subjects pooled across six
vedolizumab clinical studies: C13009 (phase 1, healthy volunteers, single
IV doses 0.2–10 mg/kg), C13002 (phase 2, UC, multiple IV 2.0/6.0/10.0
mg/kg), GEMINI 1 (phase 3, UC, 300 mg Q4W/Q8W maintenance), GEMINI 2 and
GEMINI 3 (phase 3, CD, 300 mg Q4W/Q8W maintenance), and C13004 (phase 2,
CD). Baseline demographics (Rosario 2015 Table 1 / Appendix S1 Table S1):
age 18–78 years (median 39 y in the IBD cohort), body weight 29.4–156 kg
(median 70 kg), 48% female, predominantly White (87%), 51% CD / 43% UC /
6% healthy volunteers, 50% with prior anti-TNF exposure, and 4% ADA-
positive at any time during the analysis.

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

# Source trace

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

| Equation / parameter | Value           | Source location (Rosario 2015)                                   |
|----------------------|-----------------|------------------------------------------------------------------|
| `lcl` (UC)           | `log(0.159)`    | Table 2: CLL UC = 0.159 L/day (CI 0.153–0.165)                   |
| `lcl_cd`             | `log(0.155)`    | Table 2: CLL CD = 0.155 L/day (CI 0.149–0.161)                   |
| `lvc`                | `log(3.19)`     | Table 2: Vc = 3.19 L                                             |
| `lvp`                | `log(1.65)`     | Table 2: Vp = 1.65 L                                             |
| `lq`                 | `log(0.12)`     | Table 2: Q = 0.12 L/day                                          |
| `lvmax`              | `log(0.265)`    | Table 2: Vmax = 0.265 mg/day                                     |
| `lkm`                | `log(0.964)`    | Table 2: Km = 0.964 µg/mL                                        |
| `e_wt_cl`            | `0.362`         | Table S4: weight on CLL                                          |
| `e_alb_cl`           | `-1.18`         | Table S4: albumin on CLL                                         |
| `e_calpro_cl`        | `0.0310`        | Table S4: faecal calprotectin on CLL                             |
| `e_cdai_cl`          | `-0.0515`       | Table S4: SCORE_CDAI on CLL (CD-only)                                  |
| `e_pmayo_cl`         | `0.0408`        | Table S4: partial Mayo on CLL (UC-only)                          |
| `e_age_cl`           | `-0.0346`       | Table S4: age on CLL                                             |
| `e_wt_vc`            | `0.467`         | Table S4: weight on Vc                                           |
| `e_wt_vp`            | `fixed(1)`      | Table S4: weight on Vp = 1 Fixed                                 |
| `e_wt_vmax`          | `fixed(0.75)`   | Table S4: weight on Vmax = 0.75 Fixed                            |
| `e_wt_q`             | `fixed(0.75)`   | Table S4: weight on Q = 0.75 Fixed                               |
| `e_priortnf_cl`      | `1.04`          | Table S4: prior TNF-α antagonist on CLL                          |
| `e_ada_cl`           | `1.12`          | Table S4: ADA status on CLL                                      |
| `e_conmed_aza_cl`    | `0.998`         | Table S4: AZA on CLL                                             |
| `e_conmed_mp_cl`     | `1.04`          | Table S4: MP (6-mercaptopurine) on CLL                           |
| `e_conmed_mtx_cl`    | `0.983`         | Table S4: MTX on CLL                                             |
| `e_conmed_amino_cl`  | `1.02`          | Table S4: AMINO (aminosalicylate) on CLL                         |
| `e_ibd_cd_vc`        | `1.01`          | Table S4: IBD diagnosis on Vc                                    |
| IIV `ω²_CLL`         | `0.346² = 0.1197` | Table S2/S3: %CV = 34.6 → ω = 0.346                            |
| IIV `ω²_Vc`          | `0.191² = 0.0365` | Table S2/S3: %CV = 19.1 → ω = 0.191                            |
| IIV `ω²_Vmax`        | `1.05² = 1.1025`  | Table S2/S3: %CV = 105 → ω = 1.05                              |
| `cov(CLL, Vc)`       | `+0.0374`       | Table S2: corr(CLL, Vc) = +0.566 × 0.346 × 0.191                 |
| `cov(CLL, Vmax)`     | `-0.0698`       | Table S2: corr(CLL, Vmax) = -0.192 × 0.346 × 1.05                |
| `cov(Vc, Vmax)`      | `-0.0535`       | Table S2: corr(Vc, Vmax) = -0.267 × 0.191 × 1.05                 |
| `propSd`             | `sqrt(0.0554) = 0.2354` | Table 2: σ²_prop = 0.0554 (%CV = 23.5)                  |
| Model structure      | 2-cmt parallel linear + MM | Figure 2: two-compartment model with parallel linear + Michaelis-Menten elimination from the central compartment |

Reference patient (Rosario 2015 Table 2 footnote / Figure 5 caption): 70 kg,
40 years old, albumin 4 g/dL, faecal calprotectin 700 mg/kg, SCORE_CDAI 300 (for
CD), partial Mayo 6 (for UC), UC diagnosis for Vc, no concomitant therapy,
ADA-negative, TNF-naive.

# Virtual cohort

Original observed data are not publicly available. The simulations below
use virtual populations whose covariate distributions approximate the
published phase 3 demographics (Rosario 2015 Table 1 and Appendix S1
Table S1).

```{r make-cohort}
make_cohort <- function(n,
                        n_doses = 8,
                        dosing_interval_days = 56,   # Q8W maintenance
                        obs_days_per_dose = seq(0, 56, by = 7),
                        amt_mg = 300,
                        ibd_cd_prob = 0.55,          # 55% CD, 45% UC (Table 1)
                        seed = 13243) {
  set.seed(seed)

  # Baseline covariates approximating Rosario 2015 Table 1
  WT      <- pmax(30, pmin(160, rlnorm(n, log(70), 0.22)))
  AGE     <- pmax(18, pmin(80,  rnorm(n, 39, 13)))
  ALB     <- pmax(20, pmin(55, rnorm(n, 40, 4)))        # g/dL
  SCORE_CALPRO  <- pmax(10,  exp(rnorm(n, log(700), 1.1)))         # mg/kg
  IBD_CD  <- rbinom(n, 1, ibd_cd_prob)
  # Disease-activity scores: supply SCORE_CDAI for CD, SCORE_PMAYO for UC; the gating in
  # the model sets the unused score's exponent to zero, so the supplied value
  # is irrelevant for the other cohort. We set them to the reference value to
  # make the columns self-consistent.
  SCORE_CDAI    <- ifelse(IBD_CD == 1,
                    pmax(150, pmin(600, rnorm(n, 300, 75))),
                    300)
  SCORE_PMAYO   <- ifelse(IBD_CD == 0,
                    pmax(2, pmin(9, round(rnorm(n, 6, 1.2)))),
                    6)
  PRIOR_TNF    <- rbinom(n, 1, 0.50)
  ADA_POS      <- rbinom(n, 1, 0.04)
  CONMED_AZA   <- rbinom(n, 1, 0.15)
  CONMED_MP    <- rbinom(n, 1, 0.05)
  CONMED_MTX   <- rbinom(n, 1, 0.05)
  CONMED_AMINO <- rbinom(n, 1, 0.45)

  dose_times <- seq(0, (n_doses - 1) * dosing_interval_days,
                    by = dosing_interval_days)

  pop <- data.frame(
    ID  = seq_len(n),
    WT, AGE, ALB, SCORE_CALPRO, SCORE_CDAI, SCORE_PMAYO, IBD_CD,
    PRIOR_TNF, ADA_POS, CONMED_AZA, CONMED_MP, CONMED_MTX, CONMED_AMINO
  )

  # IV infusion records: rate implements a 30-min infusion (0.0208 day).
  d_dose <- pop[rep(seq_len(n), each = length(dose_times)), ] |>
    mutate(
      TIME = rep(dose_times, times = n),
      AMT  = amt_mg,
      EVID = 1,
      CMT  = "central",
      RATE = amt_mg / (30 / (60 * 24)),
      DV   = NA_real_
    )

  # Observation records: per-dose relative sample times, across all doses
  obs_grid <- as.vector(outer(obs_days_per_dose, dose_times, "+"))
  obs_grid <- sort(unique(obs_grid))

  d_obs <- pop[rep(seq_len(n), each = length(obs_grid)), ] |>
    mutate(
      TIME = rep(obs_grid, times = n),
      AMT  = 0,
      EVID = 0,
      CMT  = "central",
      RATE = 0,
      DV   = NA_real_
    )

  bind_rows(d_dose, d_obs) |>
    arrange(ID, TIME, desc(EVID)) |>
    select(ID, TIME, AMT, EVID, CMT, RATE, DV,
           WT, AGE, ALB, SCORE_CALPRO, SCORE_CDAI, SCORE_PMAYO, IBD_CD,
           PRIOR_TNF, ADA_POS, CONMED_AZA, CONMED_MP, CONMED_MTX, CONMED_AMINO)
}
```

```{r load-model}
mod <- rxode2::rxode(readModelDb("Rosario_2015_vedolizumab"))
```

# Simulation

Phase 3 maintenance dosing: 300 mg Q8W IV × 8 doses, with weekly sampling
over the 56-week treatment period.

```{r simulate-vpc}
events_vpc <- make_cohort(n = 300)
sim_vpc <- rxode2::rxSolve(mod, events = events_vpc) |> as.data.frame()
```

## Figure 4 analogue: concentration–time VPC (phase 3 maintenance, 300 mg Q8W)

```{r fig4-vpc, fig.width=7, fig.height=4}
d_vpc <- sim_vpc |>
  group_by(time) |>
  summarise(
    Q10 = quantile(Cc, 0.10, na.rm = TRUE),
    Q50 = quantile(Cc, 0.50, na.rm = TRUE),
    Q90 = quantile(Cc, 0.90, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(d_vpc, aes(x = time, y = Q50)) +
  geom_ribbon(aes(ymin = Q10, ymax = Q90), fill = "#4682b4", alpha = 0.25) +
  geom_line(colour = "#4682b4", linewidth = 0.8) +
  geom_line(aes(y = Q10), linetype = "dashed", linewidth = 0.4) +
  geom_line(aes(y = Q90), linetype = "dashed", linewidth = 0.4) +
  scale_x_continuous(breaks = seq(0, 8 * 56, by = 56),
                     labels = seq(0, 8 * 56, by = 56) / 7) +
  labs(
    x = "Time since first dose (weeks)",
    y = expression("Vedolizumab concentration (" * mu * "g/mL)"),
    title = "Figure 4 analogue — simulated concentration VPC",
    caption = paste0(
      "Simulated 10th/50th/90th percentiles, N = 300 virtual subjects, ",
      "300 mg IV Q8W × 8 doses.\n",
      "Replicates the shape of Figure 4 (phase 3 maintenance) of Rosario 2015."
    )
  ) +
  theme_bw()
```

## Figure 5 analogue: covariate effects on CLL

Figure 5 of Rosario 2015 is a covariate forest plot showing the
multiplicative effect on steady-state linear clearance when each covariate
moves from the reference to the 5th/95th percentile (continuous) or from
0 → 1 (categorical). We compute the equivalent ratios from the model's
coefficients.

```{r fig5-forest, fig.width=7, fig.height=5}
ini <- mod$theta
get <- function(nm) unname(ini[nm])

# Continuous covariate sensitivity at 5th / 95th population percentiles.
# Percentiles approximated from the virtual cohort above.
forest_cont <- tribble(
  ~covariate,            ~low,   ~ref,   ~high,  ~exponent,
  "Weight (kg)",         49,     70,     100,    get("e_wt_cl"),
  "Albumin (g/dL)",      3.2,    4.0,    4.8,    get("e_alb_cl"),
  "Calprotectin (mg/kg)", 50,    700,    5000,   get("e_calpro_cl"),
  "Age (years)",         22,     40,     68,     get("e_age_cl"),
  "SCORE_CDAI (CD only)",      150,    300,    500,    get("e_cdai_cl"),
  "Partial Mayo (UC only)", 3,   6,      9,      get("e_pmayo_cl")
) |>
  mutate(
    ratio_low  = (low  / ref)^exponent,
    ratio_high = (high / ref)^exponent
  )

forest_cat <- tribble(
  ~covariate,                       ~ratio_low, ~ratio_high,
  "Prior TNF-α (vs naive)",         1,           get("e_priortnf_cl"),
  "ADA-positive (vs negative)",     1,           get("e_ada_cl"),
  "Concomitant AZA",                1,           get("e_conmed_aza_cl"),
  "Concomitant 6-MP",               1,           get("e_conmed_mp_cl"),
  "Concomitant MTX",                1,           get("e_conmed_mtx_cl"),
  "Concomitant aminosalicylate",    1,           get("e_conmed_amino_cl"),
  "Diagnosis CD (vs UC)",           1,           exp(get("lcl_cd") - get("lcl"))
)

forest_all <- bind_rows(
  forest_cont |> select(covariate, ratio_low, ratio_high),
  forest_cat
) |>
  mutate(covariate = factor(covariate, levels = rev(covariate)))

ggplot(forest_all, aes(y = covariate)) +
  geom_segment(aes(x = ratio_low, xend = ratio_high, yend = covariate),
               colour = "#4682b4", linewidth = 1.0) +
  geom_point(aes(x = ratio_low),  colour = "#4682b4", size = 2) +
  geom_point(aes(x = ratio_high), colour = "#4682b4", size = 2) +
  geom_vline(xintercept = 1, linetype = "dashed", colour = "grey40") +
  scale_x_log10(breaks = c(0.3, 0.5, 0.7, 1.0, 1.5, 2.0, 3.0)) +
  labs(
    x = "Multiplicative effect on linear clearance (CLL)",
    y = NULL,
    title = "Figure 5 analogue — covariate effects on CLL",
    caption = paste0(
      "Ratios computed at the 5th / 95th covariate percentiles vs the ",
      "reference.\n",
      "Reference: 70 kg, 40 y, albumin 4 g/dL, calprotectin 700 mg/kg, UC, ",
      "TNF-naive, ADA-negative, no concomitant medication."
    )
  ) +
  theme_bw()
```

# PKNCA validation

Compute NCA on the simulated typical-patient profile after a single 300 mg
IV dose and at steady state (dose 8 of a Q8W regimen). Rosario 2015 Table 4
reports the following typical exposure metrics for the 300 mg Q8W phase 3
maintenance regimen:

* Steady-state trough (week 54, pre-dose): `~9 µg/mL` for UC and CD.
* Steady-state Cmax (week 54, end-of-infusion): `~98 µg/mL` for UC and CD.
* Terminal half-life: `~25 days` (vedolizumab prescribing information
  derived from this model).

```{r pknca-single}
events_single <- make_cohort(
  n = 1, n_doses = 1,
  obs_days_per_dose = c(0.02, 0.1, 0.5, 1, 2, 4, 7, 14, 21, 28,
                        42, 56, 70, 84, 112, 140, 168)
) |>
  mutate(
    WT = 70, AGE = 40, ALB = 40, SCORE_CALPRO = 700,
    SCORE_CDAI = 300, SCORE_PMAYO = 6, IBD_CD = 0,
    PRIOR_TNF = 0, ADA_POS = 0,
    CONMED_AZA = 0, CONMED_MP = 0, CONMED_MTX = 0, CONMED_AMINO = 0
  )

mod_typical <- mod |> rxode2::zeroRe()

sim_single <- rxode2::rxSolve(mod_typical, events = events_single) |>
  as.data.frame() |>
  mutate(id = 1L, treatment = "single_300mg")

sim_nca_single <- sim_single |>
  filter(!is.na(Cc)) |>
  select(id, time, Cc, treatment)

dose_single <- events_single |>
  filter(EVID == 1) |>
  transmute(id = ID, time = TIME, amt = AMT, treatment = "single_300mg")

conc_single <- PKNCA::PKNCAconc(sim_nca_single, Cc ~ time | treatment + id)
dose_single_obj <- PKNCA::PKNCAdose(dose_single, amt ~ time | treatment + id)

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

nca_single <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(conc_single, dose_single_obj, intervals = intervals_single)
)
knitr::kable(as.data.frame(nca_single$result),
             caption = "Single-dose NCA on the typical reference-patient profile.")
```

```{r pknca-ss}
# Phase 3 maintenance steady state: 300 mg Q8W, extract the final (8th)
# dosing interval and scale time to 0 at that dose.
events_ss_base <- make_cohort(
  n = 1, n_doses = 8,
  obs_days_per_dose = c(0, 0.02, 0.1, 0.5, 1, 2, 4, 7, 14, 21, 28,
                        35, 42, 49, 56),
  dosing_interval_days = 56
) |>
  mutate(
    WT = 70, AGE = 40, ALB = 40, SCORE_CALPRO = 700,
    SCORE_CDAI = 300, SCORE_PMAYO = 6, IBD_CD = 0,
    PRIOR_TNF = 0, ADA_POS = 0,
    CONMED_AZA = 0, CONMED_MP = 0, CONMED_MTX = 0, CONMED_AMINO = 0
  )

final_dose_day <- 7 * 56   # 392 days, start of the 8th dose
extra_obs <- data.frame(
  ID = 1L,
  TIME = final_dose_day + c(70, 84, 112, 140, 168),
  AMT = 0, EVID = 0, CMT = "central", RATE = 0, DV = NA_real_,
  WT = 70, AGE = 40, ALB = 40, SCORE_CALPRO = 700,
  SCORE_CDAI = 300, SCORE_PMAYO = 6, IBD_CD = 0,
  PRIOR_TNF = 0, ADA_POS = 0,
  CONMED_AZA = 0, CONMED_MP = 0, CONMED_MTX = 0, CONMED_AMINO = 0
)
events_ss <- bind_rows(events_ss_base, extra_obs) |>
  arrange(ID, TIME, desc(EVID))

sim_ss <- rxode2::rxSolve(mod_typical, events = events_ss) |>
  as.data.frame() |>
  mutate(id = 1L, treatment = "ss_300mg_Q8W")

sim_nca_ss <- sim_ss |>
  filter(!is.na(Cc), time >= final_dose_day) |>
  mutate(time = time - final_dose_day) |>
  select(id, time, Cc, treatment)

dose_ss <- data.frame(
  id = 1L, time = 0, amt = 300, treatment = "ss_300mg_Q8W"
)

conc_ss <- PKNCA::PKNCAconc(sim_nca_ss, Cc ~ time | treatment + id)
dose_ss_obj <- PKNCA::PKNCAdose(dose_ss, amt ~ time | treatment + id)

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

nca_ss <- PKNCA::pk.nca(
  PKNCA::PKNCAdata(conc_ss, dose_ss_obj, intervals = intervals_ss)
)
knitr::kable(as.data.frame(nca_ss$result),
             caption = "Steady-state NCA on the 8th (final) dosing interval.")
```

## Comparison against published values

```{r nca-comparison}
get_param <- function(res, ppname) {
  tbl <- as.data.frame(res$result)
  val <- tbl$PPORRES[tbl$PPTESTCD == ppname]
  if (length(val) == 0) return(NA_real_)
  val[1]
}

hl_single_sim <- get_param(nca_single, "half.life")
hl_ss_sim     <- get_param(nca_ss,     "half.life")
cmax_ss_sim   <- get_param(nca_ss,     "cmax")

# Steady-state trough (C at end of dosing interval, t = 56 days post-dose 8)
c_trough_ss <- sim_ss$Cc[sim_ss$time == final_dose_day + 56]

comparison <- data.frame(
  Quantity  = c("Terminal half-life (days)",
                "Steady-state Cmax (µg/mL)",
                "Steady-state trough at week 8 (µg/mL)"),
  Published = c("~25 (vedolizumab prescribing information)",
                "~98 (Rosario 2015 Table 4 / Figure 4)",
                "~9 (Rosario 2015 Table 4 / Figure 4)"),
  Simulated = c(round(hl_ss_sim, 1),
                round(cmax_ss_sim, 1),
                round(c_trough_ss, 1))
)
knitr::kable(comparison,
             caption = "Simulated vs. published reference-patient exposures.")
```

The simulated steady-state trough, Cmax, and terminal half-life should all
fall within about 20% of the published reference-patient values. The
short-dose apparent half-life is shorter than the true terminal half-life
because the non-linear Michaelis-Menten elimination pathway contributes
meaningfully at low concentrations (Km = 0.964 µg/mL), accelerating the
terminal decline once linear-pathway clearance is no longer dominant.

# Assumptions and deviations

- **Two-typical-values switch for linear clearance.** Rosario 2015 Table 2
  reports separate typical CLL estimates for UC (0.159 L/day) and CD
  (0.155 L/day) patients. The packaged model implements this as
  `cl_typ = exp(lcl) * (1 - IBD_CD) + exp(lcl_cd) * IBD_CD` with a
  single shared IIV term `etalcl` (as the paper describes). This
  structure is the simplest faithful representation of the published
  parameterisation; an alternative ratio form
  (`exp(lcl) * (CLL_CD_ratio)^IBD_CD`) is numerically equivalent.
- **Disease-activity gating.** Partial Mayo score and SCORE_CDAI appear in the
  final model only for the diagnoses they are defined for. We gate the
  power exponents by the `IBD_CD` indicator:
  `(SCORE_PMAYO/6)^(e_pmayo_cl * (1 - IBD_CD))` and
  `(SCORE_CDAI/300)^(e_cdai_cl * IBD_CD)`. This reduces to 1 for the other
  cohort irrespective of the supplied score. For convenience, simulations
  can fill the inactive column with the reference value.
- **Time-varying ADA simplified to binary.** Rosario 2015 reports only
  4% of the pool as ADA-positive at any time and used the binary
  positivity indicator in the final model; ADA-titre sensitivity was
  not statistically significant. The packaged model therefore uses
  `ADA_POS` only; a titre effect is not present in the source.
- **IIV on Vp, Q, Km.** Rosario 2015 fixed IIV to 0 for these three
  parameters (Table S2). The packaged `ini()` block omits etas for Vp, Q,
  and Km accordingly — their typical values are shared across subjects.
- **Residual error on the un-transformed concentration scale.** Rosario
  2015 reports σ²_prop = 0.0554 as a proportional-error variance
  (Table 2). The packaged model uses `propSd = sqrt(0.0554) ≈ 0.2354`
  with `Cc ~ prop(propSd)` which encodes this in nlmixr2's standard
  proportional-error form.
- **Virtual covariate distributions.** Exact observed covariate
  distributions are not published. The demo cohort uses log-normal WT
  (median 70 kg, CV 22%), normal AGE (39 ± 13 y), normal ALB
  (4.0 ± 0.4 g/dL), log-normal SCORE_CALPRO (median 700 mg/kg, log-SD 1.1),
  55%/45% CD/UC split, and binary probabilities for the concomitant-
  medication and immunogenicity indicators approximately matching
  Rosario 2015 Table 1.

# Reference

```{r ref-footer, echo = FALSE, results = "asis"}
cat("* ", mod$reference, "\n", sep = "")
```
