---
title: "Cancer vaccine + IL-12 tumor immunotherapy in mice (Parra-Guillen 2013)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Cancer vaccine + IL-12 tumor immunotherapy in mice (Parra-Guillen 2013)}
  %\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(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
set.seed(20130419)   # day Parra-Guillen 2013 first appeared online
```

# Model and source

* Citation: `r rxode2::rxode(readModelDb("ParraGuillen_2013_cyaaE7"))$reference`
* Description: `r rxode2::rxode(readModelDb("ParraGuillen_2013_cyaaE7"))$description`
* Article: <https://doi.org/10.1208/s12248-013-9483-5>

This paper develops a semi-mechanistic K-PD tumor-growth-dynamics model on a
preclinical (mouse) CyaA-E7 cancer-vaccine dataset and then re-fits the same
structural equations to a separately published IL-12 hydrodynamic-plasmid
dataset to demonstrate the model's applicability to other immune-stimulatory
therapies. The package therefore ships **two** model files that share equations
but differ in cohort-specific parameter values:

* `modellib("ParraGuillen_2013_cyaaE7")` -- the model-build dataset (TC1 tumor
  in C57BL/6 mice, single IV dose of 50 ug CyaA-E7 vaccine on one of days 4 / 7
  / 11 / 18 / 25 / 30). All structural and IAV parameters are estimated.
* `modellib("ParraGuillen_2013_il12")` -- the applicability dataset (MC38 tumor
  in C57BL/6 mice, single hydrodynamic injection of 10 ug murine-IL-12 plasmid
  on day 23). Cell-line and immunotherapeutic-kinetics parameters (Ts0, lambda,
  k1, REG50 + IAV) are re-fit; vaccine-efficacy / regulator parameters (k3, k4,
  gamma) and the mixture probability P(1) are carried fixed from the CyaA-E7
  build.

Both files use the binary covariate `MIX_VAC_RELAPSE` (1 = relapser, 0 = cure)
to gate the SVAC degradation rate. The Bernoulli probability of relapse is
1 - P(1) = 0.156, fixed.

# Population

The CyaA-E7 model-build cohort comprised 93 female C57BL/6 mice (5 weeks old)
inoculated subcutaneously with 5 x 10^5 TC1 cells (a mouse lung-epithelial line
transformed with HPV-16 E6 / E7 and activated H-Ras). The arms were PBS day 4
(n = 19), CyaA-E7 50 ug IV day 4 (n = 15), day 7 (n = 16), day 11 (n = 14), day
18 (n = 15), day 25 (n = 14), and day 30 (n = 15). A separate 34-mouse external
validation cohort (PBS day 4 n = 23; CyaA-E7 day 25 n = 11) was held out from
the fit.

The IL-12 applicability cohort comprised 33 female C57BL/6 mice inoculated
subcutaneously with 5 x 10^5 MC38 colon-adenocarcinoma cells. Arms were PBS
control (n = 12) and IL-12 plasmid 10 ug hydrodynamic injection on day 23
(n = 21).

Tumor size was reported as the mean of two perpendicular diameters; 2 mm was
the limit of quantification (61.6 % of observations were BQL and treated as
censored using the NONMEM Laplacian M3 likelihood). Population information is
also available programmatically via
`readModelDb("ParraGuillen_2013_cyaaE7")$population` and
`readModelDb("ParraGuillen_2013_il12")$population`.

# Mathematical model

The five-ODE system (Equations 1-5 of the source, page 799):

```
dVAC/dt        = - k1 * VAC                                       (Eq. 1)
dTRAN/dt       =   k1 * VAC  - k1 * TRAN                          (Eq. 2)
dSVAC/dt       =   k1 * TRAN - k2 * SVAC                          (Eq. 3)
dTs/dt         = lambda - k3 * REG50^gamma / (REG50^gamma + REG^gamma) * Ts * SVAC   (Eq. 4)
dREG/dt        = k4 * Ts - k4 * REG                               (Eq. 5)
```

Initial conditions at tumor cell inoculation (day 0): `Ts(0) = Ts0`, the rest
zero. The vaccine virtual compartment `VAC` is reset to 1 at the time of
vaccine injection (paper Methods p. 799: "VAC0 was set to an arbitrary value
of 1 due to the absence of pharmacokinetic data"). The mixture covariate
`MIX_VAC_RELAPSE` gates `k2`:

* `MIX_VAC_RELAPSE = 0` (responder / cure, P(1) = 0.844): `k2 = 0`, so SVAC
  remains constant once formed and the tumor is eradicated permanently.
* `MIX_VAC_RELAPSE = 1` (relapser, 1 - P(1) = 0.156): `k2 = k1`, so SVAC decays
  at the same rate as the vaccine virtual drug, and the tumor regrows after
  vaccine clearance.

The shared `k1` between vaccine elimination and SVAC degradation is forced by
identifiability (Methods last paragraph p. 799).

# Source trace

| Equation / parameter | CyaA-E7 value | IL-12 value | Source location |
|-|-|-|-|
| `Ts0` (initial tumor size) | 0.324 mm | 1.16e-6 mm | Table I row `Ts0` |
| `lambda` (zero-order growth) | 0.354 mm / day | 0.335 mm / day | Table I row `lambda` |
| IAV on `lambda` | 10.1 % | 19.3 % | Table I IAV columns |
| `k1` (vaccine kinetics) | 0.0907 / day | 0.189 / day | Table I row `k1` |
| `P(1)` (responder fraction) | 0.844 FIX | 0.844 FIX | Table I row `P(1)` |
| `k2_pop1` (cure) | 0 FIX | 0 FIX | Table I row `k2_pop1` |
| `k2_pop2` (relapse) | = `k1` = 0.0907 | = `k1` = 0.189 | Table I row `k2_pop2`; Methods p. 799 |
| `k3` (vaccine efficacy) | 1.08 / day | 1.08 FIX | Table I row `k3` |
| `k4` (REG dynamics) | 0.0390 / day | 0.0390 FIX | Table I row `k4` |
| `REG50` (REG half-inhibition) | 3.18 mm | 2.08 mm | Table I row `REG50` |
| IAV on `REG50` | 33.8 % | 36.1 % | Table I IAV columns |
| `gamma` (Hill steepness) | 5.24 | 5.24 FIX | Table I row `gamma` |
| Residual SD on log(Ts) | 0.206 Log(mm) | 0.168 Log(mm) | Table I row `Residual error` |
| ODE system (Eqs. 1-5), ICs | n/a | n/a | Methods p. 799, Figure 2 |
| BQL handling (M3 likelihood) | n/a | n/a | Methods 'Model Building and Selection' p. 798 |
| `MIX_VAC_RELAPSE` covariate | Bernoulli(0.156) | Bernoulli(0.156) | Discussion p. 803; Table I `P(1)` |

# Simulation helpers

```{r helpers}
# Build an event table for one CyaA-E7 (or IL-12) treatment group. The vaccine
# enters the `vac` compartment as a unit bolus at `dose_day`; if `dose_day` is
# NA the arm is a no-dose control. `mix_relapse` is the per-mouse latent class
# (0 = responder/cure, 1 = relapser). `id_offset` keeps subject IDs disjoint
# across arms so bind_rows-ed cohorts don't collide.
make_cohort <- function(n, dose_day = NA, mix_relapse = rep(0L, n),
                        tmax = 100, id_offset = 0L) {
  stopifnot(length(mix_relapse) == n)
  obs <- tidyr::expand_grid(
    id   = id_offset + seq_len(n),
    time = seq(0, tmax, by = 1)
  )
  obs$evid <- 0L
  obs$amt  <- 0
  obs$cmt  <- "tumor_size"
  obs$MIX_VAC_RELAPSE <- rep(mix_relapse, each = tmax + 1L)
  if (!is.na(dose_day)) {
    dose <- data.frame(
      id   = id_offset + seq_len(n),
      time = dose_day,
      evid = 1L,
      amt  = 1,
      cmt  = "vac",
      MIX_VAC_RELAPSE = mix_relapse
    )
    rbind(dose, obs)
  } else {
    obs
  }
}
```

# CyaA-E7: typical responder vs. relapser individual predictions

The original paper Figure 3 plots individual model predictions for one
responder (light grey) and one relapser (dark grey) mouse per CyaA-E7 dosing
group. We reproduce the qualitative behavior with the typical-individual model
(no inter-animal variability, residual variability zeroed); responder
trajectories drop to BQL and stay there, while relapser trajectories regrow
after the vaccine signal decays.

```{r cyaaE7-typical, fig.height = 6}
mod_cyaa     <- readModelDb("ParraGuillen_2013_cyaaE7")
mod_cyaa_typ <- rxode2::zeroRe(mod_cyaa)

cyaa_days <- c(4, 7, 11, 18, 25, 30)
typ_traj <- bind_rows(lapply(cyaa_days, function(d) {
  ev_resp <- make_cohort(n = 1, dose_day = d, mix_relapse = 0L, tmax = 100)
  ev_rel  <- make_cohort(n = 1, dose_day = d, mix_relapse = 1L, tmax = 100,
                         id_offset = 1L)
  ev <- rbind(ev_resp, ev_rel)
  s <- as.data.frame(rxode2::rxSolve(mod_cyaa_typ, ev,
                                     keep = "MIX_VAC_RELAPSE"))
  s$dose_day <- d
  s$class    <- ifelse(s$MIX_VAC_RELAPSE == 1L, "relapse", "cure")
  s
}))

ggplot(typ_traj, aes(time, tumor_size, colour = class, group = interaction(id, dose_day))) +
  geom_hline(yintercept = 2, linetype = "dashed", colour = "red", alpha = 0.5) +
  geom_line(linewidth = 0.6) +
  facet_wrap(~ paste0("Vaccine day ", dose_day)) +
  scale_colour_manual(values = c("cure" = "grey60", "relapse" = "grey20")) +
  labs(x = "Time (day)", y = "Tumor size (mm)",
       title = "CyaA-E7 -- typical responder vs relapser by dosing day",
       caption = paste0(
         "Replicates Figure 3 of Parra-Guillen et al. (2013); 2 mm dashed",
         " line is the limit of quantification.")) +
  theme_bw()
```

Vaccine administration at day 4, 7, or 11 (left column) produces a sharp
drop in tumor size with the cure subpopulation reaching BQL and the
relapser bouncing back. At day 18 and beyond the regulator compartment has
accumulated enough to attenuate the vaccine-efficacy term (Eq. 4), so the
cure trajectory is shallower and the relapser barely deflects -- the source
narrative describes this as "tumour size profiles were very similar to those
obtained from the control group" at day 30 (Results, page 799).

# CyaA-E7: VPC-style population predictions

For Figure 4, the source authors simulated 1000 datasets matching the
study design and plotted the 50th percentile with 90% prediction interval.
We reproduce the same VPC envelope at one dosing day with a smaller (but
still informative) simulation, with the mixture-class assignment drawn
Bernoulli per mouse.

```{r cyaaE7-vpc, fig.height = 4.5}
sim_one_day <- function(model, dose_day, n = 200, tmax = 100,
                        p_relapse = 1 - 0.844, id_offset = 0L) {
  set.seed(dose_day + id_offset)
  mix <- as.integer(stats::runif(n) < p_relapse)
  ev  <- make_cohort(n = n, dose_day = dose_day, mix_relapse = mix,
                     tmax = tmax, id_offset = id_offset)
  out <- as.data.frame(rxode2::rxSolve(model, ev, nSub = 1L,
                                       keep = "MIX_VAC_RELAPSE"))
  out$dose_day <- dose_day
  out
}

cyaa_vpc_days <- c(4, 11, 18, 25)
cyaa_vpc <- bind_rows(lapply(seq_along(cyaa_vpc_days), function(i) {
  sim_one_day(mod_cyaa, dose_day = cyaa_vpc_days[i],
              n = 200, id_offset = i * 1000L)
}))

cyaa_vpc_summary <- cyaa_vpc |>
  group_by(dose_day, time) |>
  summarise(
    q05 = quantile(tumor_size, 0.05, na.rm = TRUE),
    q50 = quantile(tumor_size, 0.50, na.rm = TRUE),
    q95 = quantile(tumor_size, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(cyaa_vpc_summary, aes(time, q50)) +
  geom_hline(yintercept = 2, linetype = "dashed", colour = "red", alpha = 0.5) +
  geom_ribbon(aes(ymin = q05, ymax = q95), fill = "grey80", alpha = 0.6) +
  geom_line(linewidth = 0.7) +
  geom_vline(aes(xintercept = dose_day), linetype = "dotted",
             colour = "goldenrod") +
  facet_wrap(~ paste0("Vaccine day ", dose_day)) +
  labs(x = "Time (day)", y = "Tumor size (mm)",
       title = "CyaA-E7 -- simulated VPC by dosing day",
       caption = paste0(
         "Replicates Figure 4 of Parra-Guillen et al. (2013); ribbon = ",
         "5-95th percentile, line = median; gold dotted line = vaccine day.")) +
  theme_bw()
```

The widening percentile envelope across vaccine-administration arms reflects
the empirical mixture: at day 4 most mice are cured (median collapses to
BQL within ~20 days); at day 25 the regulator-driven resistance limits the
vaccine effect (median tracks the control trajectory).

# CyaA-E7: probability of cure by dosing day

Figure 5 of the source plots the simulated probability of cure (percent of
mice with predicted Ts <= 2 mm at the end of the experiment) for each
dosing group. We reproduce that summary with a larger simulation per arm.

```{r cyaaE7-cure}
end_of_study <- 60   # day used by the source for the "end-of-study" cure metric
cyaa_cure <- bind_rows(lapply(seq_along(cyaa_days), function(i) {
  sim_one_day(mod_cyaa, dose_day = cyaa_days[i], n = 200, tmax = end_of_study,
              id_offset = 10000L + i * 1000L)
})) |>
  filter(time == end_of_study) |>
  group_by(dose_day) |>
  summarise(prob_cure = mean(tumor_size <= 2, na.rm = TRUE), .groups = "drop")

cyaa_cure |>
  dplyr::rename(
    "Vaccine dose day"  = dose_day,
    "Simulated P(cure)" = prob_cure
  ) |>
  knitr::kable(
    caption = paste0("Simulated probability of cure at end of study (",
                     end_of_study, " days) -- CyaA-E7."))

ggplot(cyaa_cure, aes(dose_day, prob_cure)) +
  geom_col(width = 1.2, fill = "grey60") +
  scale_y_continuous(limits = c(0, 1)) +
  labs(x = "Vaccine dose day",
       y = "Simulated probability of cure",
       title = "CyaA-E7 -- probability of cure by dosing day",
       caption = paste0("Replicates Figure 5 of Parra-Guillen et al. ",
                        "(2013); probability of cure = fraction with ",
                        "tumor_size <= 2 mm at the end of study.")) +
  theme_bw()
```

The simulated cure rate is at or near the responder ceiling of ~84 %
(the fixed mixture probability `P(1)`) for early vaccine administration
days, declining as resistance builds for later doses -- the dominant
mechanism captured by the model.

# IL-12 applicability: typical and population predictions

The IL-12 cohort received a single hydrodynamic plasmid injection on day
23. The simulation below shows the typical responder vs relapser
trajectories and the VPC envelope; both reproduce the qualitative pattern
in Figure 4 lower-right panel of the source.

```{r il12-sim, fig.height = 4.5}
mod_il12     <- readModelDb("ParraGuillen_2013_il12")
mod_il12_typ <- rxode2::zeroRe(mod_il12)

# Typical individuals
il12_typ <- bind_rows(lapply(c(0L, 1L), function(m) {
  ev <- make_cohort(n = 1, dose_day = 23, mix_relapse = m, tmax = 60,
                    id_offset = m + 1L)
  s <- as.data.frame(rxode2::rxSolve(mod_il12_typ, ev,
                                     keep = "MIX_VAC_RELAPSE"))
  s$class <- ifelse(m == 1L, "relapse", "cure")
  s
}))

# Population VPC
il12_vpc <- sim_one_day(mod_il12, dose_day = 23, n = 200, tmax = 60,
                       id_offset = 20000L)
il12_vpc_sum <- il12_vpc |>
  group_by(time) |>
  summarise(
    q05 = quantile(tumor_size, 0.05, na.rm = TRUE),
    q50 = quantile(tumor_size, 0.50, na.rm = TRUE),
    q95 = quantile(tumor_size, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(il12_vpc_sum, aes(time, q50)) +
  geom_hline(yintercept = 2, linetype = "dashed", colour = "red", alpha = 0.5) +
  geom_ribbon(aes(ymin = q05, ymax = q95), fill = "grey80", alpha = 0.6) +
  geom_line(linewidth = 0.7) +
  geom_line(data = il12_typ, aes(time, tumor_size, colour = class, group = class),
            inherit.aes = FALSE, linewidth = 0.8) +
  geom_vline(xintercept = 23, linetype = "dotted", colour = "goldenrod") +
  scale_colour_manual(values = c("cure" = "steelblue", "relapse" = "firebrick")) +
  labs(x = "Time (day)", y = "Tumor size (mm)",
       title = "IL-12 day 23 -- typical (lines) + VPC (ribbon)",
       caption = paste0("Replicates Figure 4 lower-right of Parra-Guillen ",
                        "et al. (2013); gold dotted line = IL-12 day.")) +
  theme_bw()
```

# IL-12: probability of cure

```{r il12-cure}
il12_cure <- il12_vpc |>
  filter(time == 60) |>
  summarise(prob_cure = mean(tumor_size <= 2, na.rm = TRUE))
il12_cure |>
  dplyr::rename("Simulated P(cure) at day 60" = prob_cure) |>
  knitr::kable(
    caption = "IL-12 day 23 simulated probability of cure at day 60.")
```

The simulated cure rate sits close to the responder ceiling (84 %), as
expected for a vaccine-administration scheme that precedes substantial
resistance accumulation under the cell-line-specific REG50 = 2.08 mm.

# Assumptions and deviations

* **Mixture-class assignment.** `MIX_VAC_RELAPSE` is supplied per subject in
  the simulation event table. Typical-value figures (Figure 3 reproductions)
  set `mix_relapse` to a fixed 0 or 1. Population figures (Figure 4 / 5
  reproductions) draw `mix_relapse ~ Bernoulli(0.156)` per mouse, matching the
  fixed mixture probability `P(1) = 0.844` reported by the paper.
* **BQL handling at the simulation stage.** The model itself does not impose
  the 2 mm BQL censoring; simulated `tumor_size` values can go below 2 mm and
  even slightly negative under extreme parameter combinations. The
  probability-of-cure calculation thresholds at `tumor_size <= 2 mm` to mirror
  the paper's definition (Methods 'Visual predictive checks' / 'Numerical
  Predictive Checks' subsections, p. 800).
* **Residual error encoding.** The paper reports an additive residual on
  log-transformed tumor diameter (Methods 'Model Building and Selection',
  p. 798); the model file encodes this as `tumor_size ~ lnorm(expSd)`. The two
  forms are mathematically equivalent (additive on `log(Ts)` <=> log-normal on
  `Ts`).
* **Cohort sizes in the vignette.** The source authors simulated 1000
  populations per study design (Methods 'Model Evaluation and Validation',
  p. 800); the vignette uses 200 per arm to keep the render time within the
  pkgdown CI budget (5 minutes). The qualitative shape of the VPC envelope is
  preserved; tighter percentile bands can be obtained by raising `n` in the
  `sim_one_day()` calls above.
* **Validation dataset.** The 34-mouse external-validation cohort (PBS day 4,
  CyaA-E7 day 25) is described in the source but is not separately reproduced
  here because it shares the structural model with the training arms; the
  vaccine-day-25 panel in the CyaA-E7 VPC figure exercises the same
  conditions.
* **External datasets.** The original observed data are not redistributed
  with `nlmixr2lib`; the figures above use simulated cohorts only. The
  Berraondo 2007 (CyaA-E7) and Medina-Echeverz 2014 (IL-12) experimental data
  remain in the original publications.
