---
title: "Atezolizumab + cibisatamab colorectal cancer QSP (Anbari 2023)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Atezolizumab + cibisatamab colorectal cancer QSP (Anbari 2023)}
  %\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::rxode2(readModelDb("Anbari_2023_atezolizumab_cibisatamab_qsp"))$reference`
* Description: `r rxode2::rxode2(readModelDb("Anbari_2023_atezolizumab_cibisatamab_qsp"))$description`
* Article: <https://doi.org/10.3389/fphar.2023.1163432>
* Supplementary Material (narrative, Eqs 1-7, Figure S1): <https://www.frontiersin.org/articles/10.3389/fphar.2023.1163432/full#supplementary-material>

This is a deterministic quantitative systems pharmacology (QSP) platform, not a
population PK model. The authors built it in SimBiology and published the
complete model export as Supplementary Tables S2-S7: compartments (S2), species
and initial amounts (S3), 210 parameters (S4), 38 algebraic rules (S5), 158
reactions (S6), and the virtual-patient parameter distributions (S7). Everything
in `inst/modeldb/specificDrugs/Anbari_2023_atezolizumab_cibisatamab_qsp.R` is a
transcription of those tables; the few places where the published tables are
internally inconsistent are listed in
[Assumptions and deviations](#assumptions-and-deviations).

There are no etas and no residual-error model. The authors did not fit
inter-individual variability: they generated a 500-patient virtual cohort by
Latin hypercube sampling of the parameter distributions in Supplementary Table
S7. Encoding a fabricated eta structure would misrepresent the source, so the
model is packaged as a typical-individual mechanism model, following the same
convention as `Lu_2014_sglt_qsp`.

```{r loadmodel}
mod <- rxode2::rxode2(readModelDb("Anbari_2023_atezolizumab_cibisatamab_qsp"))
c(states = length(mod$state), parameters = length(mod$theta))
```

# Population

```{r population, echo = FALSE}
knitr::kable(
  data.frame(
    Field = c("Species", "Virtual subjects", "Calibration studies", "Disease state", "Dose range"),
    Value = c(
      "human (in silico virtual cohort)",
      "500 sampled, screened on tumour diameter, blood T cell density and Teff:Treg ratio",
      "NCT02324257 (n = 31, cibisatamab mono); NCT02650713 (n = 25, combination); IMblaze370 / NCT02788279 (atezolizumab mono)",
      "metastatic colorectal cancer",
      "cibisatamab 0-100 mg QW to Q3W; atezolizumab 1200 mg Q3W"
    )
  ),
  caption = "Virtual population (Anbari 2023, Sections 2.2 and 3.1)."
)
```

The cohort is entirely in silico. Baseline tumour diameter is uniform on
0.5-5 cm (Supplementary Table S7); the reference patient simulated throughout
this vignette uses the Table S4 baseline of 2 cm.

# Model structure

Four physiological compartments (central, peripheral, tumour, tumour-draining
lymph node), two antigen-presenting-cell compartments (endosomal volume and
endosomal / cell surface areas) and four immunological-synapse compartments
(T-cell:cancer-cell and T-cell:APC for the checkpoint module, and the two
corresponding TCE synapses).

```{r structure, echo = FALSE}
knitr::kable(
  data.frame(
    Module = c(
      "Tumour growth", "T cell dynamics", "Antigen presentation",
      "Immune checkpoint", "T cell engager", "Myeloid suppression",
      "Cytokines", "Antibody PK"
    ),
    Content = c(
      "Logistic growth of two cancer clones (C1, C2), first-order death, dead-cell clearance",
      "Naive and activated CD4+ (Treg, Th) and CD8+ (Teff) trafficking across all four compartments, thymic output, peripheral and lymph-node proliferation, exhaustion",
      "Antigen release from dying cells, APC uptake, endosomal proteolysis, peptide-MHC binding, surface translocation, TCR ligation (Hill)",
      "PD-1 / PD-L1 / PD-L2 binding in both synapses, IFN-gamma-induced PD-L1 up-regulation, bivalent atezolizumab blockade with cross-arm binding",
      "CEA-cibisatamab and CD3-cibisatamab dimers, CEA-cibisatamab-CD3 trimer formation on Teff and Treg synapses, Hill-driven enhanced killing",
      "CCL2-driven MDSC recruitment, arginase-I and nitric oxide release, T cell inhibition",
      "IL-2, TGF-beta, IFN-gamma, CCL2 secretion / degradation; Th-to-Treg transdifferentiation",
      "Four-compartment distribution with lymphatic drainage and central clearance for both antibodies (Supplementary Eqs 4-7, Table S6 R92-R97 and R145-R150)"
    )
  ),
  caption = "Modules of the Anbari 2023 QSP platform (Section 2.1 and Figure 1)."
)
```

The tumour compartment volume is not a constant: Supplementary Table S5 rule 1
makes it a repeated assignment over the cancer-cell, exhausted-T-cell and
tumour-infiltrating-T-cell counts. The packaged model therefore integrates
**amounts** for every species and derives each concentration algebraically as
amount / compartment measure, which is what SimBiology does internally and
which avoids introducing a spurious dilution term.

# Source trace

```{r sourcetrace, echo = FALSE}
knitr::kable(
  data.frame(
    Component = c(
      "Compartment capacities (V_C, V_P, V_LN, V_e, A_e, A_s, synapses)",
      "Species and initial amounts (90 states)",
      "Parameters (181 of 210 used; see Errata for the 17 unused)",
      "Algebraic and initial-assignment rules (38)",
      "Reactions / rate laws (158)",
      "TCE binding equations (narrative form)",
      "Antibody PK equations (narrative form)",
      "Reference-state values used for validation below",
      "Virtual-patient parameter distributions",
      "Simulated antibody concentration profiles",
      "kon_CD3_cibis2, kon_CEA_cibis2",
      "Atezolizumab molar mass (145 kDa)"
    ),
    `Source location` = c(
      "Supplementary Table S2",
      "Supplementary Table S3",
      "Supplementary Table S4",
      "Supplementary Table S5",
      "Supplementary Table S6",
      "Supplementary Material Eqs 1, 2-1, 2-2, 3-1, 3-2 and Table S1",
      "Supplementary Material Eqs 4-7",
      "Supplementary Table S4 (stored values of the rule outputs)",
      "Supplementary Table S7",
      "Supplementary Figure S1",
      "Ma 2020b (doi:10.1136/jitc-2020-001141) supplementary SBML export",
      "Wang 2021 (doi:10.1136/jitc-2020-002100) supplement, PK module"
    ),
    check.names = FALSE
  ),
  caption = "Provenance of every equation and parameter in the packaged model."
)
```

Each parameter line in the model file carries a trailing comment giving the
Table S4 name together with the published value **in the published units**, so
the unit conversion into the model's internal system (day, litre, nanomole,
molecule, cell) can be audited line by line.

# Initialisation: untreated tumour growth

The model is seeded with a single cancer cell (Supplementary Table S3) and grown
until it reaches the Table S4 baseline diameter of 2 cm, which is the
pre-treatment state the virtual clinical trials start from (Section 2.2).

```{r burnin}
grow <- rxode2::rxSolve(
  mod, rxode2::et(seq(0, 4600, by = 10)),
  atol = 1e-10, rtol = 1e-8, maxsteps = 1e6
) |>
  as.data.frame() |>
  mutate(diameter_cm = (6 * vol_V_T / pi)^(1 / 3) * 10)

t_start <- grow$time[which.min(abs(grow$diameter_cm - 2))]
t_start
```

```{r growthfig, echo = FALSE, fig.cap = "Untreated growth of the reference tumour from a single cell to the 2 cm pre-treatment diameter."}
grow |>
  filter(time <= t_start + 200) |>
  select(time, diameter_cm) |>
  ggplot(aes(time, diameter_cm)) +
  geom_line(linewidth = 0.8, colour = "#1f6feb") +
  geom_hline(yintercept = 2, linetype = "dashed", colour = "grey40") +
  labs(x = "Time (days)", y = "Tumour diameter (cm)") +
  theme_bw()
```

## Reference-state check against Supplementary Table S4

Supplementary Table S4 lists the rule outputs (`C_total`, `T_total`, the Hill
functions, the activated-generation numbers) alongside the true parameters.
Those entries are not inputs -- they are the values SimBiology had stored for
the calibrated reference state. They therefore act as a published, quantitative
check on the whole algebraic layer and on the immune, checkpoint and myeloid
modules simultaneously. This is the strongest validation the source supports.

```{r refstate, echo = FALSE}
ref <- c(
  C_total = 4.2136432818e8, T_total = 1.9942962767e8, Tregs_ = 2.3705305212e7,
  H_APC = 0.99982004566, H_mAPC = 0.99992788319, H_APCh = 0.99982004566,
  H_PD1_C1 = 0.99709180338, H_MDSC_C1 = 0.50840301354,
  H_TGF_reg = 0.32462605492, H_TGF_CTL = 0.19376334213,
  N_aT = 12.067318305, N_aT0 = 10.139867858, N_aTh = 10.139867858,
  H_P1 = 0.74442055179
)
row0 <- grow[which.min(abs(grow$diameter_cm - 2)), ]
cmp <- data.frame(
  Quantity   = names(ref),
  `Table S4` = signif(ref, 6),
  Model      = signif(unlist(row0[names(ref)]), 6),
  `Diff (%)` = round(100 * (unlist(row0[names(ref)]) - ref) / ref, 2),
  check.names = FALSE, row.names = NULL
)
knitr::kable(cmp, caption = "Model reference state at a 2 cm tumour versus the values stored in Supplementary Table S4.")
```

Every quantity agrees to within 5%, and the eight Hill / generation-number
terms -- which depend on the antigen-presentation, checkpoint, MDSC, TGF-beta
and IL-2 sub-models -- agree to within 2.5%. The residual differences are
consistent with the reference state being a snapshot on a slowly moving
trajectory rather than an equilibrium.

# Antibody pharmacokinetics

Supplementary Figure S1 shows the simulated concentration of each antibody in
all four compartments. Neither antibody is consumed by its binding reactions
(in Table S6 the drug appears only in the rate law, never as a reactant), so
both PK sub-systems are strictly linear and the profile expressed in mass units
is invariant to the molar mass assumed when converting the mg dose.

```{r pk}
dose_nmol <- function(mg, kda) mg * 1e-3 / (kda * 1e3) * 1e9
MW_ATEZO <- 145           # kDa, Wang 2021 supplement
MW_CIBIS <- 145           # carrier only; the ug/mL profile does not depend on it

ev_pk <- rxode2::et(seq(0, t_start + 21, by = 0.05))
ev_pk <- rxode2::et(ev_pk, amt = dose_nmol(1200, MW_ATEZO), time = t_start,
                    cmt = "q_V_C_atezo")
ev_pk <- rxode2::et(ev_pk, amt = dose_nmol(60, MW_CIBIS),
                    time = seq(t_start, t_start + 14, by = 7), cmt = "q_V_C_cibis")

pk <- rxode2::rxSolve(mod, ev_pk, atol = 1e-10, rtol = 1e-8, maxsteps = 1e7) |>
  as.data.frame() |>
  filter(time >= t_start) |>
  mutate(day = time - t_start)
```

```{r pkfig, echo = FALSE, fig.height = 6, fig.cap = "Replicates Supplementary Figure S1: atezolizumab 1200 mg Q3W (top) and cibisatamab 60 mg QW (bottom) in the central, tumour, tumour-draining lymph node and peripheral compartments."}
pk_long <- bind_rows(
  pk |> transmute(day, Drug = "Atezolizumab 1200 mg Q3W",
                  Central = x_V_C_atezo * MW_ATEZO, Tumour = x_V_T_atezo * MW_ATEZO,
                  TdLN = x_V_LN_atezo * MW_ATEZO, Peripheral = x_V_P_atezo * MW_ATEZO),
  pk |> transmute(day, Drug = "Cibisatamab 60 mg QW",
                  Central = x_V_C_cibis * MW_CIBIS, Tumour = x_V_T_cibis * MW_CIBIS,
                  TdLN = x_V_LN_cibis * MW_CIBIS, Peripheral = x_V_P_cibis * MW_CIBIS)
) |>
  pivot_longer(c(Central, Tumour, TdLN, Peripheral),
               names_to = "Compartment", values_to = "conc")

ggplot(pk_long, aes(day, conc, colour = Compartment)) +
  geom_line(linewidth = 0.7) +
  facet_wrap(~Drug, ncol = 1, scales = "free_y") +
  labs(x = "Time (days)", y = expression(Concentration~(mu*g/mL))) +
  theme_bw()
```

```{r pkcmp, echo = FALSE}
pk1 <- pk |> filter(day <= 21)
knitr::kable(
  data.frame(
    Feature = c(
      "Atezolizumab central, day 0",
      "Atezolizumab central, day 21",
      "Atezolizumab tumour : central at day 21",
      "Atezolizumab peripheral : central at day 21",
      "Cibisatamab central peak",
      "Cibisatamab central trough (day 7)"
    ),
    Model = c(
      sprintf("%.0f ug/mL", pk1$x_V_C_atezo[1] * MW_ATEZO),
      sprintf("%.0f ug/mL", tail(pk1$x_V_C_atezo, 1) * MW_ATEZO),
      sprintf("%.2f", tail(pk1$x_V_T_atezo, 1) / tail(pk1$x_V_C_atezo, 1)),
      sprintf("%.3f", tail(pk1$x_V_P_atezo, 1) / tail(pk1$x_V_C_atezo, 1)),
      sprintf("%.1f ug/mL", max(pk1$x_V_C_cibis) * MW_CIBIS),
      sprintf("%.2f ug/mL", pk1$x_V_C_cibis[which.min(abs(pk1$day - 6.9))] * MW_CIBIS)
    ),
    `Figure S1` = c("~348 ug/mL", "~73 ug/mL", "~0.41", "~0.12", "~17.6 ug/mL", "~0.3 ug/mL"),
    check.names = FALSE
  ),
  caption = "Model versus values digitised from Supplementary Figure S1."
)
```

The shape, the decay kinetics and the tumour-to-central and
peripheral-to-central ratios reproduce Figure S1 well. The absolute levels are
lower than the figure by a roughly constant factor, and the lymph-node
compartment sits further below the figure than the other three. That gap is a
property of the source, not of the transcription: the narrative PK equations
(Supplementary Eqs 4-7) contain no interstitial-volume-fraction terms at all,
while the Table S6 rate laws divide every driving concentration by `gamma`, and
neither reading reproduces the absolute levels in Figure S1. The packaged model
follows Table S6, which is the actual SimBiology export. See
[Assumptions and deviations](#assumptions-and-deviations).

## Non-compartmental analysis of the simulated profiles

```{r nca}
conc_df <- pk |>
  filter(day <= 21) |>
  transmute(id = 1L, treatment = "Atezolizumab 1200 mg",
            time = day, Cc = x_V_C_atezo * MW_ATEZO) |>
  filter(!is.na(Cc))

dose_df <- data.frame(id = 1L, treatment = "Atezolizumab 1200 mg",
                      time = 0, amt = 1200, duration = 0)

o_conc <- PKNCA::PKNCAconc(conc_df, Cc ~ time | treatment / id)
o_dose <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id,
                           duration = "duration")
res_nca <- PKNCA::pk.nca(PKNCA::PKNCAdata(o_conc, o_dose))
```

```{r ncatab, echo = FALSE}
nca_tab <- as.data.frame(res_nca) |>
  filter(PPTESTCD %in% c("cmax", "tmax", "auclast", "half.life", "cl.last")) |>
  select(PPTESTCD, PPORRES) |>
  mutate(PPORRES = signif(PPORRES, 4)) |>
  rename("NCA parameter" = PPTESTCD, "Simulated value" = PPORRES)
knitr::kable(nca_tab, caption = "NCA of the simulated atezolizumab central profile over one 21-day cycle.")
```

The model has no published NCA table to compare against, so this serves as an
internal consistency check on the PK sub-model: the terminal half-life should be
governed by the slow peripheral redistribution rather than by
`CL_atezo / V_C = 0.324 / 5 = 0.065` per day (a 10.7-day half-life), and the
simulated value is consistent with that.

# Treatment comparison

The paper's headline result (Section 3.2 and Figure 6) is that concurrent
combination therapy gives a better response than either monotherapy. The
reference patient below is the typical individual, not the 500-patient cohort,
so it reproduces the ordering rather than the overall response rates in Table 1
(which are cohort quantities requiring the Table S7 Latin hypercube sample).

```{r arms}
HORIZON <- 400
Da <- dose_nmol(1200, MW_ATEZO)
Dc <- dose_nmol(60, MW_CIBIS)
arms <- list(
  "No treatment"    = list(),
  "Cibisatamab QW"  = list(list(cmt = "q_V_C_cibis", amt = Dc, ii = 7)),
  "Atezolizumab Q3W" = list(list(cmt = "q_V_C_atezo", amt = Da, ii = 21)),
  "Combination"     = list(list(cmt = "q_V_C_cibis", amt = Dc, ii = 7),
                           list(cmt = "q_V_C_atezo", amt = Da, ii = 21))
)
tx <- lapply(names(arms), function(nm) {
  ev <- rxode2::et(seq(0, t_start + HORIZON, by = 5))
  for (d in arms[[nm]]) {
    ev <- rxode2::et(ev, amt = d$amt,
                     time = seq(t_start, t_start + HORIZON, by = d$ii), cmt = d$cmt)
  }
  rxode2::rxSolve(mod, ev, atol = 1e-10, rtol = 1e-8, maxsteps = 1e7) |>
    as.data.frame() |>
    filter(time >= t_start) |>
    mutate(day = time - t_start, arm = nm,
           diameter_cm = (6 * vol_V_T / pi)^(1 / 3) * 10)
}) |>
  bind_rows()
```

```{r armfig, echo = FALSE, fig.cap = "Replicates the treatment ordering of Figure 3B / Figure 6: tumour diameter over 400 days of monotherapy and concurrent combination therapy in the reference patient."}
ggplot(tx, aes(day, diameter_cm, colour = arm)) +
  geom_line(linewidth = 0.8) +
  labs(x = "Days since treatment start", y = "Tumour diameter (cm)", colour = NULL) +
  theme_bw()
```

```{r armtab, echo = FALSE}
tx |>
  group_by(arm) |>
  summarise(
    `Diameter at week 8 (cm)` = round(diameter_cm[which.min(abs(day - 56))], 3),
    `Diameter at day 400 (cm)` = round(dplyr::last(diameter_cm), 3),
    `Best change from baseline (%)` = round(100 * (min(diameter_cm) - dplyr::first(diameter_cm)) /
                                              dplyr::first(diameter_cm), 1),
    .groups = "drop"
  ) |>
  arrange(`Diameter at day 400 (cm)`) |>
  knitr::kable(caption = "Reference-patient response by arm.")
```

Concurrent combination therapy gives the smallest tumour at both the week-8
first-follow-up landmark and at the end of treatment, and is the only arm that
produces any regression from baseline in the reference patient. Both
monotherapies slow growth relative to no treatment. This is the ordering
reported in Section 3.2.

# Assumptions and deviations

Everything below is a documented departure from a literal reading of the
published tables, or a value the published tables do not supply. No parameter
value was taken from outside the on-disk sources.

1. **Table S6 omits the neo-antigen MHC translocation reaction.** The
   self-antigen (P0) and neo-antigen (P1) presentation modules are exact mirrors
   of each other (R73/R83, R74/R84, R75/R85, R76/R86, R77/R87, R78/R88, R79/R89,
   R80/R90) with one exception: R81, `A_e.M1p0 -> A_s.M1p0`, has no P1
   counterpart. Without it the surface complex `A_s.M1p1` has a sink (R90) and
   no source, so it decays to zero, `H_P1` collapses to about 1e-9, and no CD8 T
   cell response can occur. That contradicts Supplementary Table S4, which
   stores `pTCR_p1_MHC_tot = 7.7055e-5` and `H_P1 = 0.74442` for the reference
   state. The mirror reaction `A_e.M1p1 -> A_s.M1p1` with rate
   `kout * A_e.M1p1 * A_e` is therefore added, and is labelled `v81b` in the
   model. With it, `H_P1` reproduces to 2.5%.

2. **Table S6 reactions 72 and 82 drop the antigen-load factor on two terms.**
   In the antigen-deposition rate laws the first killing term carries the
   `P0_C1` (respectively `P1_C1`) factor that converts a cell-death rate into an
   antigen-release rate, but the `k_C_BTcell` and `k_C_BTreg` terms do not. As
   printed the expression adds mol/day to cell/day, which is dimensionally
   impossible. The missing `P0_C1` / `P0_C2` (`P1_C1` / `P1_C2`) factors are
   restored by analogy with the adjacent term and with R61 / R62.

3. **`kon_CD3_cibis2` and `kon_CEA_cibis2` are not in Table S4.** Table S6
   reactions 152, 154, 156 and 158 use these second-step two-dimensional
   on-rates, but Table S4 lists only `kon_CD3_TCE` and `kon_CEA_TCE`. The
   upstream SimBiology export of Ma 2020b, which this paper states the TCE
   binding module is based on, gives `kon_CD3_cibis2 = 3333.33` and
   `kon_CEA_cibis2 = 333.33` in 1/(molarity*nanometer*second) against
   `kon_CD3_cibis = 10000` and `kon_CEA_cibis = 1000` in 1/(molarity*second) --
   that is, `kon_X_cibis2 = kon_X_cibis / D_syn` with the synapse gap
   `D_syn = 3 nm`, which is also this paper's Table S4 value. The same relation
   is applied to this paper's on-rates, giving `10000 / 3` and `20000 / 3`.

4. **Six Table S6 symbols are renames of Table S4 entries.** `SA_cell` is
   `SA_Ccell`; `k_cl_atezo` and `k_cl_cibis` are `CL_atezo` and `CL_cibis`
   (confirmed by Supplementary Eq 4, where the clearance term is `CL * A_C`);
   and `kon_CEA_cibis`, `koff_CEA_cibis`, `kon_CD3_cibis`, `koff_CD3_cibis` are
   the `*_TCE` entries. Each is flagged in the model file.

5. **Avogadro conversions are explicit.** Twelve surface-surface bimolecular
   reactions (102, 103, 105, 107, 112, 113, 115, 117, 152, 154, 156, 158) and
   the four endosomal peptide-MHC reactions (78, 79, 88, 89) mix a molar
   association constant with molecule-per-area surface densities. SimBiology
   inserts the molecule-to-mole conversion silently; the packaged model divides
   by `NAVG_nmol` explicitly at each of those points. This is exactly the
   `k_on / (D_syn * N)` grouping written out in Supplementary Eqs 1 and 3.

6. **Cross-compartment transport is scaled by the tumour volume.** Five rate
   laws (95, 96, 146, 148, 149) are concentration rates whose reactant and
   product sit in different compartments. Scaling each species by its own
   compartment would not conserve mass. Supplementary Eqs 4-7 write every
   lymphatic-drainage term as `q_LD * V_T * A`, for both the losing and the
   gaining compartment, so the tumour volume is used for all five. The same
   applies to `q_T_cibis`, which Table S4 gives in 1/second where its
   atezolizumab counterpart `q_T_atezo` is in mL/second; multiplying by the
   tumour volume makes it the volumetric flow that Eq 6 requires.

7. **Figure S1 absolute levels are not reproduced.** As described above, the
   narrative PK equations and the Table S6 rate laws disagree about the
   interstitial volume fractions, and neither reproduces the absolute
   concentrations in Figure S1. The packaged model follows Table S6. Relative
   behaviour (decay kinetics, tumour-to-central and peripheral-to-central
   ratios) matches; the lymph-node level is the poorest agreement.

8. **Cibisatamab molar mass is not reported anywhere on disk.** The main paper,
   the supplement, and the upstream Ma 2020a, Ma 2020b, Wang 2021 and
   Jafarnejad 2019 supplements all give cibisatamab doses in mg without a molar
   mass, and the two primary sources for the molecule (Van de Vyver 2021,
   Bacac 2016) are not open access. The model therefore doses in nmol, its
   native unit, and no molar mass is encoded. The mg-to-nmol conversion used in
   this vignette assumes 145 kDa purely as a carrier: the mass-unit PK profile
   is invariant to that choice because the PK is linear, but the efficacy
   simulation in the treatment-comparison section is not. Read the
   treatment-comparison figure as the ordering of the arms, which is what the
   paper reports, rather than as a calibrated dose-response.

9. **Seventeen Table S4 parameters are unused.** `A_syn`, `T_syn`, `SA_APC`,
   `SA_Ccell`, `CL_atezo`, `CL_cibis`, `kon_CEA_TCE`, `koff_CEA_TCE`,
   `kon_CD3_TCE`, `koff_CD3_TCE` and `D_syn` are duplicated under the names the
   rate laws actually use (items 3 and 4 above) or are the compartment
   capacities of Table S2. `N_avg` is folded into `NAVG_nmol`. `H_PD1_APC`,
   `N_costim`, `T_CD28_total` and `T_CTLA4_syn` are inherited from the parent
   models and are referenced by no rule or reaction in Table S5 or S6.
   `initial_tumour_diameter` is a simulation-protocol parameter, not a model
   parameter: it is the diameter at which treatment begins, and is used above to
   pick the treatment start time rather than inside `model()`.

10. **No inter-individual variability or residual error.** The source reports
    none; see the note under [Model and source](#model-and-source). To build the
    virtual cohort, sample the 21 parameters of Supplementary Table S7 and pass
    them to `rxode2::rxSolve()` as a parameter data frame.

11. **`covariateData` is empty.** The model takes no subject-level covariate
    columns. Every input is either a fixed parameter or a dosing event; the
    between-patient variation in the source is variation in model parameters,
    not in observed covariates.

12. **State names deliberately depart from the canonical compartment register.**
    `checkModelConventions()` reports 90 warnings, one per state, and no errors.
    Each state is a (compartment, species) pair in the SimBiology sense, and is
    named `q_<compartment>_<species>` directly from the Parent and Name columns
    of Supplementary Table S3 -- for example `q_V_T_C1` is species `C1` in
    compartment `V_T`, and `q_syn_T_C1_PD1_PDL1` is the PD-1:PD-L1 complex in
    the T-cell:cancer-cell synapse. Mapping these onto canonical names is not
    possible (there are four drug states per antibody, so no single state can be
    `central`), and registering 90 model-specific names as new canonicals would
    not help any other model. The systematic naming keeps every state
    one-to-one auditable against Table S3, which is the property that matters
    for a transcription of this size. The `q_` prefix marks the state as an
    amount; the derived `x_<compartment>_<species>` variables in `model()` are
    the corresponding SimBiology-valued concentrations or surface densities.
