---
title: "Pacmilimab conditionally activated anti-PD-L1 TNBC QSP (Ippolito 2024)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Pacmilimab conditionally activated anti-PD-L1 TNBC QSP (Ippolito 2024)}
  %\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

```{r model}
# built once and reused throughout: this is a 240-state model and each build
# costs about a minute
mod <- rxode2::rxode2(readModelDb("Ippolito_2024_pacmilimab_qsp"))
length(mod$state)
```

* Citation: `r mod$reference`
* Article: <https://doi.org/10.1002/psp4.13060>
* Appendix S1 (Eqs S1-S68) and Tables S1-S2: supporting information of the article

Pacmilimab is a Probody therapeutic (PbTx): an anti-PD-L1 antibody whose two
active sites are each covered by a mask held on by a protease-cleavable linker.
The mask shifts reversibly to reveal the active site (equilibrium constant
`K_M`), and is removed irreversibly by tumour-associated proteases (rate
`k_cvg`), which are far more active inside the tumour than outside it. The
paper embeds this dynamic in the Popel-laboratory QSP-IO platform previously
validated for atezolizumab in triple-negative breast cancer, and adds an
out-of-synapse module so PD-L1 outside the immunological synapse can also bind
drug.

This is a **deterministic mechanism model**, not a population model. The
authors generate virtual patients by Latin hypercube sampling of the 29
parameter distributions in Supplementary Table S1; there are no etas and no
residual error model, so everything below is a typical-value simulation.

# Population

```{r population, echo = FALSE}
pop <- readModelDb("Ippolito_2024_pacmilimab_qsp")
pop <- eval(body(pop)[[which(vapply(as.list(body(pop)),
  function(e) is.call(e) && identical(e[[2]], as.name("population")), TRUE))]][[3]])
knitr::kable(data.frame(Field = names(pop), Value = unlist(lapply(pop, as.character))),
             row.names = FALSE)
```

The clinical anchors are the pacmilimab phase I study (16 patients with TNBC,
objective response rate 7%) and the single-dose PK of Stroh 2021 at 0.1, 0.3,
1, 3 and 10 mg/kg. The regimen selected in that trial, and simulated
throughout the paper, is 10 mg/kg every two weeks for 400 days.

# Model structure

| Layer | Content |
|---|---|
| Physiological compartments | central, peripheral, tumour, tumour-draining lymph node |
| Antigen-presentation compartments | APC endosome (`V_e`), APC surface (`A_s`, `A_e`) |
| Synapse compartments | T cell-cancer cell, T cell-APC, macrophage-cancer cell |
| Out-of-synapse compartments (new in this paper) | cancer cell, T cell, macrophage, APC |
| Cell populations | two cancer clones, naive/activated/effector CD8 and CD4 T cells, Tregs, exhausted T cells, APC and mature APC, MDSC, M1 and M2 macrophages |
| Soluble mediators | IL-2, IL-10, IL-12, IFN-gamma, TGF-beta, CCL2, nitric oxide, arginase-I, angiogenic factor |
| Checkpoints | PD-1/PD-L1/PD-L2, CTLA-4/CD28/CD80/CD86, CD47/SIRP-alpha |
| Probody | six free states per compartment, nine bound states at each of nine PD-L1 bearing sites |

States are integrated as **amounts** (`q_<compartment>_<species>`) and the
declared-unit value (concentration, surface density or cell count) is derived
algebraically as `x_<compartment>_<species>`. This is what SimBiology does
internally and it removes the dilution question for the tumour compartment,
whose volume `vol_V_T` is a state-dependent repeated assignment over the cell
counts it contains.

The six free Probody states are named for the status of the two active sites:
`m` masked, `o` open (mask present but revealing the site), `c` cleaved. So
`P_mm` is fully masked (the form injected as pacmilimab), `P_cc` is fully
cleaved (equivalent to the naked anti-PD-L1 antibody), and `P_co` has one site
cleaved and one open. Per Figure 4a of the paper, the **active** states are
every state except `P_mm`.

# Source trace

| Model element | Source |
|---|---|
| Probody state transitions (unmasking, cleavage), 4 compartments | Ippolito Appendix S1 Eqs S51-S56 |
| Probody transport and central clearance | Ippolito Appendix S1 Eqs S30-S33 |
| Probody mono- and bivalent PD-L1 binding, 9 states per site | Ippolito Appendix S1 Eqs S57-S65 |
| Out-of-synapse PD-L1 binding and turnover, 4 compartments | Ippolito Appendix S1 Eqs S66-S68 |
| Unmasking equilibrium `K_M` = 1/175, unmasking rate `k_o` = 0.0116 /s | Ippolito Table S2 |
| Cleavage rate `k_cvg` = 1e-07 /s (tumour) | Ippolito Table S1, "Probody cleavage rate in TNBC" |
| Extra-tumoral cleavage = 1% of the tumour rate | Ippolito Appendix S1, "Model implementation" |
| All 237 remaining parameter values | Ippolito Table S2 |
| Cancer, T cell, APC, antigen, checkpoint, MDSC and macrophage ODEs | Wang 2022 iScience Tables S3-S6 (the SimBiology export of the platform this paper inherits, cited as ref. 3 of Appendix S1) |
| `kon_CD80_CD80`, `koff_CD80_CD80` | Wang 2022 Table S2 (absent from Ippolito Table S2; see Errata) |
| Antibody molar mass 145 kDa | Wang 2022 supplementary methods (atezolizumab; the paper assumes the PbTx shares its PK) |

Every parameter carries an inline comment in the model file naming its source
table and its value in the source's own units, alongside the converted value
used by the model.

```{r consts}
MW <- 145000                                   # g/mol, 145 kDa
WT <- 70                                       # kg (not reported; see Errata)
nmol_of <- function(mgkg) mgkg * WT / MW * 1e6
diam_cm <- function(v_L) (6 * v_L * 1000 / pi)^(1/3)
# rtol 1e-7 reproduces rtol 1e-8 to 7-8 significant figures on every quantity
# used below, at ~200x less cost on the stiff early growth phase
SOLVE <- list(atol = 1e-10, rtol = 1e-7, maxsteps = 1e7)
solv <- function(...) do.call(rxode2::rxSolve,
  c(list(mod, ...), SOLVE, list(returnType = "data.frame")))

FREE   <- c("P_mm", "P_mo", "P_oo", "P_cm", "P_co", "P_cc")
ACTIVE <- setdiff(FREE, "P_mm")                # Figure 4a
UNCLEAVED <- c("P_mm", "P_mo", "P_oo")         # Figure 2b assay
```

# Initialisation: untreated tumour growth

The paper's virtual-patient protocol simulates each patient from a small
cancer-cell population until the tumour reaches the pre-treatment size, then
starts dosing. Table S1 gives a median pre-treatment diameter of 2.5 cm.

```{r grow}
grow <- solv(rxode2::et(seq(0, 320, by = 5)))
grow$diam <- diam_cm(grow$vol_V_T)
day_start <- grow$time[which.min(abs(grow$diam - 2.5))]
day_start
```

```{r growfig, echo = FALSE, fig.cap = "Untreated Gompertzian growth from the Table S3 initial cancer-cell population to the 2.5 cm median pre-treatment diameter of Table S1."}
ggplot(grow, aes(time, diam)) +
  geom_line(linewidth = 0.8) +
  geom_hline(yintercept = 2.5, linetype = 2, colour = "grey40") +
  geom_vline(xintercept = day_start, linetype = 3, colour = "grey40") +
  labs(x = "Day", y = "Tumour diameter (cm)") +
  theme_bw()
```

The state at that point is the pre-treatment patient used for every arm below.

```{r ref}
ref <- grow[which.min(abs(grow$diam - 2.5)), ]
inits <- unlist(ref[mod$state]); names(inits) <- mod$state
inits[inits < 0] <- 0
```

```{r reftab, echo = FALSE}
hn <- c("H_PD1_C1", "H_PD1_APC", "H_CD28_APC", "H_APC", "H_P0", "H_P1", "H_TGFb",
        "H_MDSC", "H_IL10", "H_IL12", "H_SIRPa", "H_PD1_M", "H_Mac_C", "N_aT", "N_aT0")
knitr::kable(data.frame(Quantity = c(hn, "Teff/Treg", "M1/M2"),
  Value = signif(c(unlist(ref[hn]), ref$x_V_T_T1 / ref$Tregs_,
                   ref$x_V_T_Mac_M1 / ref$x_V_T_Mac_M2), 5)),
  row.names = FALSE, caption = "Pre-treatment reference state of the typical patient.")
```

Unlike the sibling Anbari 2023 export, Ippolito's Table S2 stores SimBiology's
*default* values for the repeated-assignment targets (`H_APC` 0.5, `H_PD1_C1`
0.9, `N_aT` 1, `C_total` 0) rather than converged reference-state outputs, so
this table is a description of the model's own pre-treatment state rather than
a check against published numbers. Those rows are model **outputs** and are
deliberately not carried into `ini()`; see Errata.

# Probody mass conservation

Two conservation identities follow from Eqs S51-S65 and hold independently of
any parameter value, so they are the sharpest available test that the Probody
equations were transcribed correctly.

**Free states.** Unmasking and cleavage only move drug between the six free
states and transport only moves it between compartments, so free drug plus
cumulative central clearance must equal the dose at all times. Note that free
drug is deliberately *not* decremented by surface binding: Appendix S1 assumes
"an abundance of probody and antibody in each compartment compared to the ones
interacting with the cell surface", exactly as the parent model's own aPDL1
reactions do. Free-plus-bound is therefore not a conserved quantity.

```{r masscons}
mk_ev <- function(target, amt, grid, ii = NULL, addl = NULL) {
  d <- if (is.null(ii)) rxode2::et(amt = amt, time = 0, cmt = target)
       else rxode2::et(amt = amt, time = 0, cmt = target, ii = ii, addl = addl)
  rxode2::add.sampling(d, grid)
}
dose10 <- nmol_of(10)
pk_grid <- seq(0, 84, by = 0.25)
pk_pb <- solv(mk_ev("q_V_C_P_mm", dose10, pk_grid), inits = inits)

free_drug <- function(s) rowSums(s[, as.vector(outer(c("V_C", "V_P", "V_T", "V_LN"),
  FREE, function(a, b) paste0("q_", a, "_", b))), drop = FALSE])
clear_rate <- function(s) rowSums(s[, paste0("vpb_cl_", FREE), drop = FALSE])
cum_trap <- function(t, y) c(0, cumsum(diff(t) * (head(y, -1) + tail(y, -1)) / 2))

bal <- free_drug(pk_pb) + cum_trap(pk_pb$time, clear_rate(pk_pb))
max(abs(bal / dose10 - 1))
stopifnot(max(abs(bal / dose10 - 1)) < 1e-3)
```

**Bound states.** With association switched off, the cleavage and unmasking
terms of Eqs S57-S65 must cancel exactly, leaving dissociation of the
monovalent complexes as the only sink. (A bivalent complex that loses one bond
becomes a monovalent complex, so its `2*koff` term is a transfer within the
bound pool, not a sink from it.)

```{r boundcons}
BOUND <- c("P_mo_fb", "P_oo_fb", "P_cm_bf", "P_co_bf", "P_co_fb", "P_cc_fb",
           "P_oo_bb", "P_co_bb", "P_cc_bb")
site <- "syn_T_C1_PDL1"
i2 <- inits
for (b in BOUND) i2[[paste0("q_", site, "_", b)]] <- 1
bs <- solv(rxode2::et(seq(0, 0.4, by = 0.001)), inits = i2,
           params = c(kon_PDL1_aPDL1 = 0))
mono <- rowSums(bs[, paste0("q_", site, "_", BOUND[1:6]), drop = FALSE])
biv  <- rowSums(bs[, paste0("q_", site, "_", BOUND[7:9]), drop = FALSE])
tot  <- mono + biv
k    <- which(tot > 1e-6); k <- k[k > 1 & k < length(tot)]
dtot <- (tot[k + 1] - tot[k - 1]) / (bs$time[k + 1] - bs$time[k - 1])
pred <- (-bs$p_koff_PDL1_aPDL1[1] * mono)[k]
max(abs((dtot - pred) / pred))
stopifnot(max(abs((dtot - pred) / pred)) < 1e-3)
```

Both identities hold to solver tolerance. The second one is what caught three
defects in the printed equations (see Errata) - it fails on five states if
Eqs S61 and S62 are implemented exactly as printed.

# Antibody pharmacokinetics

Figure 2b of the paper calibrates the model against single-dose blood
concentrations of **uncleaved** pacmilimab at 0.1, 0.3, 1, 3 and 10 mg/kg.
Uncleaved means no site has been cleaved, i.e. `P_mm + P_mo + P_oo`.

```{r pkdose}
doses <- c(0.1, 1, 10)
pk <- lapply(doses, function(d) {
  s <- solv(mk_ev("q_V_C_P_mm", nmol_of(d), pk_grid), inits = inits)
  data.frame(id = 1L, treatment = paste0(d, " mg/kg"), dose_mgkg = d,
             time = s$time,
             Cc = rowSums(s[, paste0("q_V_C_", UNCLEAVED), drop = FALSE]) / s$m_V_C)
}) |> bind_rows()
```

```{r pkfig, echo = FALSE, fig.cap = "Replicates Figure 2b of Ippolito 2024: simulated blood uncleaved pacmilimab after a single dose."}
ggplot(pk, aes(time, Cc, colour = treatment)) +
  geom_line(linewidth = 0.8) +
  scale_y_log10() +
  labs(x = "Day", y = "Blood uncleaved pacmilimab (nmol/L)", colour = "Dose") +
  theme_bw()
```

## Non-compartmental analysis of the simulated profiles

```{r nca}
conc_df <- pk |> dplyr::filter(!is.na(Cc))
dose_df <- pk |>
  dplyr::distinct(id, treatment, dose_mgkg) |>
  dplyr::mutate(time = 0, amt = nmol_of(dose_mgkg))
o_conc <- PKNCA::PKNCAconc(conc_df, Cc ~ time | treatment / id)
o_dose <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id,
                           route = "intravascular", duration = 0)
# one explicit interval spanning the whole profile, so every parameter comes
# from the same interval (PKNCA's default puts auclast on 0-24 h and the rest
# on 0-Inf, which would mix intervals in the table below)
ivl <- data.frame(start = 0, end = max(pk_grid), cmax = TRUE, tmax = TRUE,
                  auclast = TRUE, aucinf.obs = TRUE, half.life = TRUE, cl.obs = TRUE)
res_nca <- PKNCA::pk.nca(PKNCA::PKNCAdata(o_conc, o_dose, intervals = ivl))
```

```{r ncatab, echo = FALSE}
nca_raw <- as.data.frame(res_nca) |>
  dplyr::filter(PPTESTCD %in% c("cmax", "tmax", "auclast", "aucinf.obs",
                                "half.life", "cl.obs")) |>
  dplyr::select(treatment, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |>
  dplyr::left_join(dplyr::distinct(pk, treatment, dose_mgkg), by = "treatment")
knitr::kable(dplyr::mutate(nca_raw, dplyr::across(where(is.numeric), ~ signif(.x, 4))),
             row.names = FALSE,
             caption = "PKNCA summary of the simulated uncleaved-pacmilimab profiles.")
```

The paper publishes no NCA table for pacmilimab, so there are no reference
values to compare against. What *is* checkable is that the free-Probody
sub-system is strictly linear - drug appears in the transport, unmasking and
cleavage terms only as a first-order reactant, and surface binding does not
consume it - so every exposure metric must be exactly dose-proportional.

```{r linearity}
spread <- function(x) max(x) / min(x) - 1
lin <- nca_raw |>
  dplyr::mutate(cmax_per_mgkg = cmax / dose_mgkg,
                auclast_per_mgkg = auclast / dose_mgkg,
                aucinf_per_mgkg = aucinf.obs / dose_mgkg)
c(cmax = spread(lin$cmax_per_mgkg),
  auclast = spread(lin$auclast_per_mgkg),
  aucinf = spread(lin$aucinf_per_mgkg),
  half.life = spread(lin$half.life),
  cl.obs = spread(lin$cl.obs))
stopifnot(
  spread(lin$cmax_per_mgkg) < 1e-8,
  spread(lin$auclast_per_mgkg) < 1e-5,
  spread(lin$aucinf_per_mgkg) < 1e-5,
  spread(lin$half.life) < 1e-5,
  spread(lin$cl.obs) < 1e-5
)
```

Dose-normalised `Cmax`, `AUClast` and `AUCinf`, the terminal half-life and the
clearance are identical across a 100-fold dose range, as they must be.

Only a small fraction of the dose is cleaved systemically, which is the direct
consequence of the paper's assumption that extra-tumoral cleavage runs at 1% of
the tumour rate:

```{r cleaved}
cleaved_frac <- tail(rowSums(pk_pb[, paste0("q_V_C_", c("P_cm", "P_co", "P_cc")),
                                   drop = FALSE]), 1) /
  tail(rowSums(pk_pb[, paste0("q_V_C_", FREE), drop = FALSE]), 1)
signif(cleaved_frac, 3)
```

# Treatment comparison

The paper simulates three arms over 400 days: untreated, unmasked anti-PD-L1
monotherapy, and pacmilimab monotherapy, all at 10 mg/kg every two weeks. The
unmasked arm is produced by dosing the fully cleaved state `P_cc` - Appendix S1
defines it as "injecting a fully cleaved (i.e., without both masks) antibody" -
so a single model covers all three arms.

```{r arms}
tr_grid <- seq(0, 400, by = 5)
arm <- function(target) solv(mk_ev(target, dose10, tr_grid, ii = 14, addl = 28),
                             inits = inits)
untreated <- solv(rxode2::et(tr_grid), inits = inits)
unmasked  <- arm("q_V_C_P_cc")
pbtx      <- arm("q_V_C_P_mm")

arms <- bind_rows(
  data.frame(arm = "untreated",  time = untreated$time, vol = untreated$vol_V_T,
             H_PD1 = untreated$H_PD1_C1),
  data.frame(arm = "unmasked",   time = unmasked$time,  vol = unmasked$vol_V_T,
             H_PD1 = unmasked$H_PD1_C1),
  data.frame(arm = "pacmilimab", time = pbtx$time,      vol = pbtx$vol_V_T,
             H_PD1 = pbtx$H_PD1_C1)) |>
  group_by(arm) |>
  mutate(diam = diam_cm(vol), pct_change = 100 * (vol / first(vol) - 1)) |>
  ungroup()
```

```{r armfig, echo = FALSE, fig.cap = "Replicates the tumour-size comparison of Figure 3 of Ippolito 2024 for the typical patient: pacmilimab is intermediate between the untreated and unmasked arms."}
ggplot(arms, aes(time, diam, colour = arm)) +
  geom_line(linewidth = 0.8) +
  labs(x = "Day", y = "Tumour diameter (cm)", colour = NULL) +
  theme_bw()
```

```{r armtab, echo = FALSE}
arms |>
  group_by(arm) |>
  summarise(`Diameter day 0 (cm)` = signif(first(diam), 4),
            `Diameter day 400 (cm)` = signif(last(diam), 4),
            `Volume change (%)` = signif(last(pct_change), 4), .groups = "drop") |>
  dplyr::rename(Arm = arm) |>
  knitr::kable(row.names = FALSE)
```

```{r armorder}
final <- arms |> group_by(arm) |> summarise(v = last(vol), .groups = "drop")
stopifnot(
  final$v[final$arm == "unmasked"] < final$v[final$arm == "pacmilimab"],
  final$v[final$arm == "pacmilimab"] < final$v[final$arm == "untreated"]
)
```

The ordering reproduces the paper's central efficacy conclusion: masking costs
some efficacy, but both antibodies are far better than no treatment. Figure S1
of the paper shows the corresponding fall in the PD-1 checkpoint Hill term, and
the same ordering appears here.

```{r hpd1, echo = FALSE}
arms |>
  group_by(arm) |>
  summarise(`H_PD1 at day 400` = signif(last(H_PD1), 5), .groups = "drop") |>
  dplyr::rename(Arm = arm) |>
  knitr::kable(row.names = FALSE)
```

# Tumour selectivity of activation

Figure 4 of the paper compares the concentration of *active* antibody (any
state with at least one open or cleaved site) between compartments, and reports
that masking raises the tumour-to-periphery and tumour-to-blood ratios by
roughly 10% while dropping the absolute peripheral exposure by more than an
order of magnitude.

```{r selectivity}
act_conc <- function(s, cp) {
  cap <- if (cp == "V_T") s$vol_V_T else s[[paste0("m_", cp)]]
  rowSums(s[, paste0("q_", cp, "_", ACTIVE), drop = FALSE]) / cap
}
act_frac <- function(s, cp) rowSums(s[, paste0("q_", cp, "_", ACTIVE), drop = FALSE]) /
                            rowSums(s[, paste0("q_", cp, "_", FREE), drop = FALSE])

sel <- data.frame(
  Arm = c("unmasked", "pacmilimab"),
  `Active in tumour (nmol/L)`     = c(tail(act_conc(unmasked, "V_T"), 1),
                                      tail(act_conc(pbtx, "V_T"), 1)),
  `Active in periphery (nmol/L)`  = c(tail(act_conc(unmasked, "V_P"), 1),
                                      tail(act_conc(pbtx, "V_P"), 1)),
  `Active fraction, tumour`       = c(tail(act_frac(unmasked, "V_T"), 1),
                                      tail(act_frac(pbtx, "V_T"), 1)),
  `Active fraction, periphery`    = c(tail(act_frac(unmasked, "V_P"), 1),
                                      tail(act_frac(pbtx, "V_P"), 1)),
  `Active fraction, blood`        = c(tail(act_frac(unmasked, "V_C"), 1),
                                      tail(act_frac(pbtx, "V_C"), 1)),
  check.names = FALSE)
sel$`Selectivity T/P` <- sel$`Active fraction, tumour` / sel$`Active fraction, periphery`
sel$`Selectivity T/blood` <- sel$`Active fraction, tumour` / sel$`Active fraction, blood`
```

```{r seltab, echo = FALSE}
knitr::kable(dplyr::mutate(sel, dplyr::across(where(is.numeric), ~ signif(.x, 4))),
             row.names = FALSE)
```

```{r selassert}
periph_fold <- sel$`Active in periphery (nmol/L)`[1] / sel$`Active in periphery (nmol/L)`[2]
periph_fold
stopifnot(
  periph_fold > 10,                       # paper: "more than an order of magnitude"
  sel$`Selectivity T/P`[2] > 1.10,        # masking concentrates activation in the TME
  sel$`Selectivity T/blood`[2] > 1.10,
  abs(sel$`Selectivity T/P`[1] - 1) < 1e-9  # the unmasked antibody is active everywhere
)
```

The active fraction of circulating pacmilimab is about 2%, and it is roughly
22% higher inside the tumour than in the periphery or the blood - the
conditional-activation effect the molecule is designed for. The unmasked
antibody is by construction 100% active in every compartment, so its
selectivity is exactly 1. Absolute active exposure in the periphery is
`r signif(periph_fold, 3)`-fold lower for pacmilimab, comfortably reproducing
the paper's "more than an order of magnitude" claim.

Note that the *raw concentration* ratio `[A]_T/[A]_P` plotted in Figures 4b and
4d is not reproduced in the same direction at typical-value parameters, because
the two arms end the 400 days with different tumour volumes and the tumour
concentration is diluted accordingly. The paper averages that ratio over 10,000
virtual patients, where the between-arm difference in tumour size is much
smaller. The volume-free active-fraction comparison above isolates the
mechanism the paper is describing and does reproduce it.

# Assumptions and deviations

* **Antibody molar mass and body weight.** Doses are reported in mg/kg. The
  molar mass is not given for pacmilimab; the model assumes the 145 kDa of
  atezolizumab reported in the Wang 2022 supplementary methods, consistent with
  the paper's own statement that the PbTx shares atezolizumab's PK. Body weight
  is not reported anywhere; 70 kg is used. Both only set the scale of the dose
  in nmol: the free-Probody sub-system is exactly linear, so every ratio,
  fraction and dose-normalised metric in this vignette is invariant to them.
* **Three defects in the printed Probody equations.** Requiring that cleavage
  and dissociation conserve antibody across Eqs S57-S65 determines all three
  uniquely:
  1. Eq S61 loses `k_cvg [P_co,bf]`; it must be `k_cvg [P_co,fb]`, otherwise
     `P_co,fb` has no cleavage sink and `P_co,bf` has two.
  2. Eq S61 gains `koff [P_oo,bb]`; it must be `koff [P_co,bb]`, since
     `P_oo,bb` already returns `2*koff` to `P_oo,fb` in Eq S58 while the
     `2*koff` sink of `P_co,bb` in Eq S64 is otherwise unmatched.
  3. Eq S62 carries a dangling `+ k_cvg` with no species and writes its own
     dissociation sink as `koff [P_co,fb]`. The dangling term is
     `k_cvg [P_cm,bf]` - the only cleavage flux in Eq S59 with no destination -
     and the sink is `koff [P_cc,fb]`.
  The bound-state conservation check above fails on five states without these
  corrections.
* **Rule outputs are not parameters.** SimBiology stores the last computed
  value of every repeated-assignment rule alongside the true parameters, so
  Table S2 lists `C_total`, `T_total`, `M_total`, `C_max`, every `H_*` Hill
  term, the `N_a*` generation numbers, `R_Tcell`, `Tregs_` and the `pTCR_*`
  totals as though they were constants. They are outputs and are computed by
  their rules; carrying them into `ini()` would pin `C_total` to 0 and freeze
  the tumour volume. Their stored values here are SimBiology defaults, not
  converged reference-state outputs, so they provide no independent check.
* **The nab-paclitaxel arm is excluded.** Ippolito's Table S2 omits all twelve
  nab-paclitaxel PK/PD parameters of the parent model and pins
  `k_C1_therapy = k_C2_therapy = 0`. The four `NabP` species, the three rules
  that reference them and the seven chemotherapy reactions are therefore
  dropped, together with the `k_C_resist` clone-switching reaction (resistance
  in the parent model is resistance to nab-paclitaxel). Cancer clone `C2`
  consequently stays at its initial value of 0.
* **`kon_CD80_CD80` and `koff_CD80_CD80`** are used by the inherited
  CD80/CD28/CTLA-4 network but are absent from Ippolito's Table S2. They are
  taken from the on-disk upstream Wang 2022 Table S2, the structural source
  this paper cites for the whole checkpoint module. No value was invented.
* **State count.** The paper reports 218 ODEs and 45 algebraic rules; this
  implementation has 240 ODEs. The difference is in how widely Eqs S57-S65 are
  instantiated. Appendix S1 states that they are "implemented for both the
  synapse and out-of-synapse compartments" and that the out-of-synapse module
  declares one compartment for each of the four PD-L1 expressing cell types, so
  the implementation here instantiates the nine bound states at all nine PD-L1
  bearing sites, which is the reading the equations support. A narrower reading
  would reach 218 but the paper does not say which sites to omit.
* **Free drug is not consumed by surface binding.** This is the paper's
  explicit abundance assumption and matches the parent model's own aPDL1
  reactions. Free-plus-bound antibody is therefore not conserved.
* **No IIV and no residual error.** The paper generates virtual patients by
  Latin hypercube sampling of the Table S1 distributions rather than by
  estimating a covariance matrix, so no etas and no error model are encoded.
  Table S1 is reproduced in the model file comments for anyone wanting to build
  the cohort.
* **Cohort size.** The paper simulates 10,000 virtual patients; this vignette
  runs typical-value arms only, which is what the deterministic structure
  supports within a vignette time budget.
* **Non-canonical state names.** The `q_<compartment>_<species>` state names
  are machine-derived from the source tables so that every state is one-to-one
  auditable against Wang 2022 Table S3 and Appendix S1. `checkModelConventions()`
  reports them as non-canonical compartment names; registering 240
  model-specific canonical names would not help any other model.
