---
title: "Etrolizumab (Moein 2022)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Etrolizumab (Moein 2022)}
  %\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("Moein_2022_etrolizumab"))
cat("* **Citation:** ", mod$reference, "\n", sep = "")
cat("* **Description:** ", mod$description, "\n", sep = "")
```

* **Article:** <https://doi.org/10.1002/psp4.12846>

# Population

The model was developed on 1263 subjects (non-missing covariates) pooled from
five clinical studies of etrolizumab in adults with moderately-to-severely
active ulcerative colitis: ABS4262g (phase I), EUCALYPTUS (phase II), HIBISCUS
I/II (phase III), HICKORY (phase III), and LAUREL (phase III). GARDENIA
(phase III) was held out for external validation. Baseline demographics
(Moein 2022 Table 3): age range 18-79 years (median 38 years), weight range
38.0-216 kg (median 72.2 kg), 42% female, 45% with prior anti-TNF exposure,
23% ADA-positive. Disease extension: 53% left-sided colitis, 42%
extensive/pancolitis, 2% other, 3% missing. Baseline medians: albumin
41 g/L, CRP 4.31 mg/L, fecal calprotectin 1500 ug/g. Dosing ranged from
single doses of 0.3-10 mg/kg IV (phase I), 0.5-3 mg/kg Q4W SC (phase I),
315-420 mg loading SC (phase II), 4 mg/kg Q4W IV (phase I), to the
phase III regimen of 105 mg Q4W SC.

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

# Source trace

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

| Equation / parameter | Value            | Source location (Moein 2022)                         |
|----------------------|------------------|------------------------------------------------------|
| `lka`                | `log(0.193)`     | Table 4: ka = 0.193 /day                             |
| `lcl`                | `log(0.260)`     | Table 4: CL (at TSFD = 0) = 0.260 L/day              |
| `lvc`                | `log(2.61)`      | Table 4: Vc = 2.61 L                                 |
| `lvp`                | `log(1.77)`      | Table 4: Vp = 1.77 L                                 |
| `lq`                 | `log(0.449)`     | Table 4: Q = 0.449 L/day                             |
| `logitfdepot`        | `logit(0.712)`   | Table 4: F = 0.712                                   |
| `logitmaxred`        | `logit(0.263)`   | Table 4: Maxred = 0.263                              |
| `lonset`             | `log(4.81)`      | Table 4: Onset = 4.81 weeks                          |
| `e_wt_cl_q`          | `0.872`          | Table 4 note a: WT exponent on CL and Q              |
| `e_wt_vc_vp`         | `0.788`          | Table 4 note b: WT exponent on Vc and Vp             |
| `e_alb_cl`           | `-0.0314`        | Table 4: Albumin on CL                               |
| `e_crp_cl`           | `0.00458`        | Table 4: CRP on CL                                   |
| `e_adat_cl`          | `0.0365`         | Table 4: ADAT on CL                                  |
| `e_priortnf_cl`      | `0.0490`         | Table 4: Prior TNF on CL                             |
| `e_extpan_cl`        | `0.0816`         | Table 4: Extensive/pancolitis on CL                  |
| `e_othext_cl`        | `0.181`          | Table 4: Other disease extension on CL               |
| IIV CL (`etalcl`)    | `log(1+0.243^2)` | Table 4: IIV CL CV = 0.243                           |
| IIV Vc (`etalvc`)    | `log(1+0.252^2)` | Table 4: IIV Vc CV = 0.252                           |
| IIV Vp (`etalvp`)    | `log(1+0.262^2)` | Table 4: IIV Vp CV = 0.262                           |
| IIV F                | `0.733^2`        | Table 4: IIV F logit-SD = 0.733                      |
| IIV Maxred           | `0.597^2`        | Table 4: IIV Maxred logit-SD = 0.597                 |
| `propSd`             | `0.196`          | Table 4: Proportional residual CV = 0.196            |
| `addSd`              | `0.427`          | Table 4: Additive residual SD = 0.427 ug/mL          |
| CL time-decay        | Equation 1       | Eq. 1: CL = CL0 * (1 - Maxred * (1 - exp(-log(2)/(Onset*7) * (TSFD - TAD)))) |
| Covariate form       | `exp(theta*(Cov - Cov_ref))` | Table 4 note d (continuous)              |
| Covariate form       | `1 + theta * indicator`      | Table 4 note g (categorical)             |

Reference covariate values for the reference patient (Figure 1 caption):
72 kg WT, 41 g/L albumin, 4.23 mg/L CRP, phase III, no prior anti-TNF,
left-sided colitis, ADA-negative.

# Virtual cohort

Original observed data are not publicly available. The simulations below
use virtual populations whose covariate distributions approximate the
published phase III demographics.

```{r make-cohort}
make_cohort <- function(n, n_doses = 6, dosing_interval_days = 28,
                        obs_days_per_dose = seq(0, 28, by = 2),
                        amt_mg = 105,
                        seed = 12846) {
  set.seed(seed)

  # Baseline covariates approximating the phase III population
  WT  <- pmax(35, pmin(200, rlnorm(n, log(72), 0.22)))
  ALB <- pmax(25, pmin(55, rnorm(n, 41, 4)))
  CRP <- pmax(0.1, exp(rnorm(n, log(4.31), 1.2)))
  PRIOR_TNF <- rbinom(n, 1, 0.45)
  disext <- sample(c("left", "extpan", "other"), n,
                   replace = TRUE, prob = c(0.55, 0.43, 0.02))
  DISEXT_EP    <- as.integer(disext == "extpan")
  DISEXT_OTHER <- as.integer(disext == "other")
  # Paper assumes ADA-negative until week 4; use a time-invariant titer of 0
  # for the primary cohort and a 97.5th-percentile value for a sensitivity arm
  ADA_TITER <- rep(0, n)

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

  pop <- data.frame(
    ID  = seq_len(n),
    WT, ALB, CRP, ADA_TITER, PRIOR_TNF, DISEXT_EP, DISEXT_OTHER
  )

  # Dose records (depot, SC injections Q4W)
  d_dose <- pop[rep(seq_len(n), each = length(dose_times)), ] |>
    mutate(
      TIME = rep(dose_times, times = n),
      AMT  = amt_mg,
      EVID = 1,
      CMT  = "depot",
      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",
      DV   = NA_real_
    )

  bind_rows(d_dose, d_obs) |>
    arrange(ID, TIME, desc(EVID)) |>
    select(ID, TIME, AMT, EVID, CMT, DV,
           WT, ALB, CRP, ADA_TITER, PRIOR_TNF, DISEXT_EP, DISEXT_OTHER)
}
```

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

# Simulation

Phase III dosing: 105 mg Q4W SC x 6 doses, with observations every 2 days over
the 24-week treatment period.

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

# Replicate published figures

## Figure 2 analogue: concentration-time VPC (phase III dose)

```{r fig2-vpc-plot, 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, 168, by = 28),
                     labels = seq(0, 168, by = 28) / 7) +
  labs(
    x = "Time since first dose (weeks)",
    y = expression("Etrolizumab concentration (" * mu * "g/mL)"),
    title = "Figure 2 analogue - simulated concentration VPC",
    caption = paste0(
      "Simulated 10th/50th/90th percentiles, N = 300 virtual subjects, ",
      "105 mg Q4W SC x 6 doses.\n",
      "Replicates the shape of Figure 2 of Moein 2022 (phase III panels)."
    )
  ) +
  theme_bw()
```

## Figure S3 analogue: typical CL vs. time since first dose

Time-dependent CL per Equation 1 reaches its asymptote over the first
4-8 weeks, stepping down after each dose.

```{r figS3-cl-typical, fig.width=6.5, fig.height=3.5}
# Population-typical CL trajectory for the reference patient (72 kg, 41 g/L
# ALB, 4.23 mg/L CRP, no prior TNF, left-sided colitis, ADA-negative).
# zeroRe() removes IIV to get typical values.
mod_typical <- mod |> rxode2::zeroRe()
events_typ <- make_cohort(n = 1) |>
  mutate(WT = 72, ALB = 41, CRP = 4.23,
         ADA_TITER = 0, PRIOR_TNF = 0, DISEXT_EP = 0, DISEXT_OTHER = 0)
sim_typ <- rxode2::rxSolve(mod_typical, events = events_typ,
                           keep = c("WT", "ALB", "CRP")) |>
  as.data.frame()

ggplot(sim_typ, aes(x = time, y = cl)) +
  geom_line(linewidth = 0.8, colour = "#4682b4") +
  scale_x_continuous(breaks = seq(0, 168, by = 28),
                     labels = seq(0, 168, by = 28) / 7) +
  labs(
    x = "Time since first dose (weeks)",
    y = "Typical CL (L/day)",
    title = "Figure S3 analogue - typical CL vs time since first dose",
    caption = paste0(
      "Reference patient (72 kg, 41 g/L ALB, 4.23 mg/L CRP, no prior TNF, ",
      "left-sided colitis, ADA-).\n",
      "Maxred = 26% asymptotic reduction; Onset = 4.8 weeks half-life."
    )
  ) +
  theme_bw()
```

# PKNCA validation

Compute NCA on the simulated typical-value profile for the reference patient:

1. **Single-dose phase** using the first dosing interval (0-28 days) as the
   analog for the phase I single-dose pK profile.
2. **Steady-state phase** using the final dosing interval (days 140-168) where
   the time-dependent CL has reached its asymptote.

The paper reports typical terminal half-lives derived from the population PK
parameters:

* After a single dose: 13.0 days (95% CI 12.2-13.9).
* At steady state: 17.1 days (95% CI 16.1-18.3).

For a clean terminal-phase characterization we simulate a single 105 mg dose
with dense post-peak sampling, and a 6-dose steady-state regimen with one
dense post-final-dose sampling interval.

## Single-dose NCA

```{r pknca-single}
# Single 105 mg SC dose, dense sampling out to 84 days (3 half-lives past
# the expected 13-day t1/2).
events_single <- make_cohort(
  n = 1, n_doses = 1,
  obs_days_per_dose = c(0, 0.25, 0.5, 1, 2, 3, 5, 7, 10, 14, 21,
                        28, 35, 42, 49, 56, 63, 70, 77, 84)
) |>
  mutate(WT = 72, ALB = 41, CRP = 4.23,
         ADA_TITER = 0, PRIOR_TNF = 0, DISEXT_EP = 0, DISEXT_OTHER = 0)

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

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_105mg")

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-patient profile.")
```

## Steady-state NCA

```{r pknca-ss}
# Six 105 mg SC doses Q4W, with dense post-peak sampling in the final
# interval for terminal-phase characterization.
events_ss <- make_cohort(
  n = 1, n_doses = 6,
  obs_days_per_dose = c(0, 0.25, 0.5, 1, 2, 3, 5, 7, 10, 14, 21, 28)
) |>
  mutate(WT = 72, ALB = 41, CRP = 4.23,
         ADA_TITER = 0, PRIOR_TNF = 0, DISEXT_EP = 0, DISEXT_OTHER = 0)

# Add dense post-final-dose terminal-phase samples out to day 252 (84 days
# past the final dose, for t1/2 estimation at steady state)
final_dose_day <- 5 * 28
extra_obs <- data.frame(
  ID = 1L, TIME = final_dose_day + c(35, 42, 49, 56, 63, 70, 77, 84),
  AMT = 0, EVID = 0, CMT = "central", DV = NA_real_,
  WT = 72, ALB = 41, CRP = 4.23,
  ADA_TITER = 0, PRIOR_TNF = 0, DISEXT_EP = 0, DISEXT_OTHER = 0
)

events_ss <- bind_rows(events_ss, extra_obs) |>
  arrange(ID, TIME, desc(EVID))

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

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 = 105, treatment = "ss_105mg_Q4W"
)

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 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")

comparison <- data.frame(
  Quantity     = c("Terminal half-life after single dose (days)",
                   "Terminal half-life at steady state (days)"),
  Published    = c("13.0 (95% CI 12.2-13.9)",
                   "17.1 (95% CI 16.1-18.3)"),
  Simulated    = c(round(hl_single_sim, 2), round(hl_ss_sim, 2))
)
knitr::kable(comparison,
             caption = "Simulated vs. published terminal half-lives.")
```

The simulated terminal half-lives should fall within about 20% of the
published central estimates. The single-dose value is shorter than the
steady-state value because CL is at its maximum (0.260 L/day) after the
first dose and drops to approximately 0.260 * (1 - 0.263) = 0.192 L/day
at steady state, lengthening t1/2 proportionally.

# Assumptions and deviations

- **PHASE12 covariate on residual error not implemented.** Moein 2022
  Table 4 reports a -19.2% relative change in the combined residual error
  for phase I/II studies vs. phase III. nlmixr2's `add()` / `prop()`
  residual-error terms do not accept per-observation expressions, so the
  packaged model uses the phase III residual-error magnitude
  (`propSd = 0.196`, `addSd = 0.427 ug/mL`) unconditionally. For typical
  phase III simulations (the primary usage pattern) this is exact; for
  phase I/II simulations the residual error will be slightly over-stated.
  If a future release of nlmixr2 supports per-observation error
  magnitudes, this covariate can be added.
- **Time-varying ADA titer simplified.** Moein 2022 treats ADA titer as a
  time-varying covariate. The reference patient used in the source
  figures (Figure 1) is ADA-negative at least until week 4, so
  `ADA_TITER = 0` is the appropriate reference for reproducing the
  published figures. Users simulating ADA-positive patients should
  supply a time-series of titer values per subject.
- **Reference CRP.** The model centers CRP at 4.23 mg/L (the reference-
  patient value reported in the Figure 1 caption). The population
  baseline median CRP is 4.31 mg/L (Table 3); both values produce
  effectively identical exposures.
- **Virtual covariate distributions.** Exact baseline distributions are
  not published. The demo cohort uses log-normal WT (median 72 kg, CV
  22%), normal ALB (41 +/- 4 g/L), log-normal CRP (median 4.31 mg/L,
  log-SD 1.2), binary prior anti-TNF (45% probability), and a three-
  level disease extension (55% left / 43% extpan / 2% other) to
  approximate Table 3 of Moein 2022.
- **SC-only dosing.** The packaged model uses SC dosing into `depot`.
  Moein 2022 also pooled phase I IV dose-levels; for IV simulations
  the event table should dose into `central` directly and
  bioavailability `f(depot) = 1` is irrelevant. The time-dependent CL
  remains correct because `time - tad(depot)` still resolves to the
  time of the most recent dose; when dosing IV into `central`, use
  `tad(central)` or pre-compute the TSFD/TAD difference in the event
  table.

# Reference

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