Wang_2024_atezolizumab_carboplatin_nabpaclitaxel_qsp <- function() {
  description <- paste(
    "QSP. Quantitative systems pharmacology model of non-small-cell lung cancer",
    "(NSCLC) describing the triple combination of the anti-PD-L1 checkpoint",
    "inhibitor atezolizumab with carboplatin and nab-paclitaxel, extended from",
    "the published immuno-oncology QSP platform. Four physiological compartments",
    "(central, peripheral, tumour, tumour-draining lymph node) plus APC",
    "endosomal and surface compartments, two immunological-synapse compartments",
    "and a separate three-compartment Michaelis-Menten PK model for",
    "nab-paclitaxel. Covers naive and activated CD4+/CD8+ T cell trafficking and",
    "proliferation, logistic tumour growth against a dynamic vasculature-driven",
    "carrying capacity with a drug-sensitive and a chemo-resistant cancer clone,",
    "self-antigen and neo-antigen release, uptake, endosomal processing and",
    "MHC-I presentation, TCR-mediated activation, PD-1/PD-L1/PD-L2 checkpoint",
    "binding with atezolizumab blockade, CTLA-4/CD28/CD80/CD86 co-signalling,",
    "MDSC recruitment with arginase-I and nitric-oxide suppression,",
    "Th-to-Treg transdifferentiation, TGF-beta / IFN-gamma / IL-2 / CCL2",
    "dynamics, and the cytotoxic plus anti-angiogenic effects of nab-paclitaxel",
    "and carboplatin. 140 ODE states, 55 algebraic (repeated- and",
    "initial-assignment) rules and 247 parameters, translated from the",
    "SimBiology export in Supplementary Tables S2-S7. Deterministic mechanism",
    "model: the authors generated virtual patients by Latin hypercube sampling",
    "of the 26 parameter distributions in Supplementary Table S1 rather than by",
    "fitting IIV or residual error, so no etas and no error model are encoded.",
    sep = " "
  )
  reference <- paste(
    "Wang CY, Dai HR, Tan YP, Yang DH, Niu XM, Han L, Wang W, Ma LL, Julku A,",
    "Jiao Z. Development and Evaluation of a Quantitative Systems Pharmacology",
    "Model for Mechanism Interpretation and Efficacy Prediction of Atezolizumab",
    "in Combination with Carboplatin and Nab-Paclitaxel in Patients with",
    "Non-Small-Cell Lung Cancer. Pharmaceuticals (Basel). 2024;17(2):238.",
    "doi:10.3390/ph17020238.",
    sep = " "
  )
  vignette <- "Wang_2024_atezolizumab_carboplatin_nabpaclitaxel_qsp"

  units <- list(
    time          = "day",
    dosing        = "nmol",
    concentration = "nmol/L"
  )

  covariateData <- list()

  population <- list(
    species       = "human (in silico virtual cohort)",
    n_subjects    = 1000L,
    n_studies     = 1L,
    disease_state = "advanced squamous / non-squamous non-small-cell lung cancer, first line",
    dose_range    = paste(
      "Atezolizumab 1200 mg IV on day 1 of each 21-day cycle (with 840 mg Q2W",
      "and 1680 mg Q4W evaluated as alternative regimens), carboplatin at an",
      "AUC of 6 mg/mL/min IV on day 1, and nab-paclitaxel 100 mg/m2 IV on days",
      "1, 8 and 15.",
      sep = " "
    ),
    notes = paste(
      "Calibrated and evaluated against the IMpower131 trial (NCT02367794).",
      "Virtual cohorts were generated by Latin hypercube sampling of the 26",
      "parameter distributions in Supplementary Table S1 (500 patients for the",
      "objective-response-rate and duration-of-response evaluation, 1000 per",
      "regimen for the alternative-dosing comparison and the biomarker",
      "analysis). Pre-treatment tumour diameter is log-normal with a geometric",
      "mean of 2.5 cm. Treatment duration 400 days with tumour assessment every",
      "8 weeks, scored by RECIST v1.1. Immune-system parameters were assumed",
      "shared with the authors' earlier triple-negative-breast-cancer platform;",
      "only the cancer-cell dynamics parameters were re-estimated for NSCLC.",
      sep = " "
    )
  )

  ini({
    # ---- Compartment capacities (Supplementary Table S2) ----
    vol_V_C <- fixed(5); label("Capacity of compartment V_C [dm^3]")  # Table S2 (V_C = 5 liter)
    vol_V_P <- fixed(60); label("Capacity of compartment V_P [dm^3]")  # Table S2 (V_P = 60 liter)
    # V_T is a repeated assignment (Table S6 rule 1), not a constant: stored value 1e-06 milliliter
    vol_V_LN <- fixed(0.00137444678595); label("Capacity of compartment V_LN [dm^3]")  # Table S2 (V_LN = 1374.45 millimeter^3)
    vol_V_e <- fixed(4e-16); label("Capacity of compartment V_e [dm^3]")  # Table S2 (V_e = 4e-16 litre)
    vol_A_e <- fixed(1.5e-09); label("Capacity of compartment A_e [dm^2]")  # Table S2 (A_e = 15 micrometer^2)
    vol_A_s <- fixed(9e-08); label("Capacity of compartment A_s [dm^2]")  # Table S2 (A_s = 900 micrometer^2)
    vol_syn_T_C1 <- fixed(3.78e-09); label("Capacity of compartment syn_T_C1 [dm^2]")  # Table S2 (syn_T_C1 = 37.8 micrometer^2)
    vol_syn_T_APC <- fixed(3.78e-09); label("Capacity of compartment syn_T_APC [dm^2]")  # Table S2 (syn_T_APC = 37.8 micrometer^2)
    vol_V_1 <- fixed(15.8); label("Capacity of compartment V_1 [dm^3]")  # Table S2 (V_1 = 15.8 liter)
    vol_V_2 <- fixed(1650); label("Capacity of compartment V_2 [dm^3]")  # Table S2 (V_2 = 1650 liter)
    vol_V_3 <- fixed(75.4); label("Capacity of compartment V_3 [dm^3]")  # Table S2 (V_3 = 75.4 liter)

    NAVG_nmol <- fixed(6.0221409e14); label("Avogadro constant, molecule per nanomole [molecule/nmol]")  # physical constant; carries the molecule <-> nmol conversion

    # ---- Model parameters (Supplementary Table S4) ----
    # Values are converted to the canonical unit system (day / L / dm^2 / nmol / cell / ug);
    # each trailing comment records the value and unit exactly as tabulated in Table S4.
    k_cell_clear <- fixed(0.1); label("Rate of dead cell clearance from tumour compartment (estimated) [1/day]")  # Table S4 (k_cell_clear = 0.1 1/day)
    cell <- fixed(1); label("unit parameter for calculation [cell]")  # Table S4 (cell = 1 cell)
    day <- fixed(1); label("unit parameter for calculation [day]")  # Table S4 (day = 1 day)
    vol_cell <- fixed(2.52731128245e-12); label("Average volume of cancer cell Volume of a cancer cell calculated based on cancer cell diameter (Jafarnejad 2019 PMID: 31236847) [dm^3/cell]")  # Table S4 (vol_cell = 2527.31 micrometer^3/cell)
    vol_Tcell <- fixed(1.750157098e-13); label("Average volume of T cells Volume of a T cell calculated based on the average T cell diameter (Chapman 1981, PMID: 6975780) [dm^3/cell]")  # Table S4 (vol_Tcell = 175.016 micrometer^3/cell)
    Ve_T <- fixed(0.37); label("Void fraction of the tumor (Finley 2012, PMID: 22547351) [1]")  # Table S4 (Ve_T = 0.37 dimensionless)
    k_C1_growth <- fixed(0.012); label("Cancer cell growth rate (doi: 10.1208/s12248-019-0302-5) [1/day]")  # Table S4 (k_C1_growth = 0.012 1/day)
    k_C1_death <- fixed(0.0001); label("Cancer cell death rate from innate immune cells (PMID: 3156189) [1/day]")  # Table S4 (k_C1_death = 0.0001 1/day)
    initial_tumour_diameter <- fixed(0.25); label("Pre-treatment tumor diameter [dm]")  # Table S4 (initial_tumour_diameter = 2.5 centimeter)
    div_T0 <- fixed(1160000); label("T Cell Diversity (Robins 2009, PMID: 19706884) [1]")  # Table S4 (div_T0 = 1.16e+06 dimensionless)
    n_T0_clones <- fixed(250); label("Number of T cell clones (Jafarnejad 2019 PMID: 31236847) [1]")  # Table S4 (n_T0_clones = 250 dimensionless)
    q_nT0_LN_in <- fixed(0.1); label("Rate of naive T cell transport into the LN calculated based on the steady-state naive T cell level in healthy individuals (Autissier 2010, PMID: 20099249) [1/day]")  # Table S4 (q_nT0_LN_in = 0.1 1/day)
    q_T0_LN_out <- fixed(24); label("Rate of activated T cell transport out of the LN (De Boer and Perelson 1995, PMID: 7475092) [1/day]")  # Table S4 (q_T0_LN_out = 24 1/day)
    q_nT0_LN_out <- fixed(2.88); label("Rate of naive T cell transport out of the LN (Mandl 2012, PMID: 23071319) [1/day]")  # Table S4 (q_nT0_LN_out = 2.88 1/day)
    k_T0_act <- fixed(5); label("T0 activation rate (De Boer and Perelson 1995, PMID: 7475092; Robertson-Tessi 2012, PMID: 22051568) [1/day]")  # Table S4 (k_T0_act = 5 1/day)
    k_T0_pro <- fixed(1); label("T0 proliferation rate (Marchingo 2014, PMID: 25430770) [1/day]")  # Table S4 (k_T0_pro = 1 1/day)
    k_T0_death <- fixed(0.01); label("T0 death rate (De Boer and Perelson 1995, PMID: 7475092) [1/day]")  # Table S4 (k_T0_death = 0.01 1/day)
    q_T0_P_in <- fixed(3507.84); label("rate of T0 transport into the peripheral compartment calculated based on T cell transmigration rate (Zhu 1996, PMID: 8706023) and T cell adhesion density (Zhu 1996, PMID: 8706023) [1/day]")  # Table S4 (q_T0_P_in = 2.436 1/minute)
    q_T0_P_out <- fixed(24); label("rate of T0 transport out of the peripheral compartment (De Boer and Perelson 1995, PMID: 7475092) [1/day]")  # Table S4 (q_T0_P_out = 24 1/day)
    q_T0_T_in <- fixed(83.52); label("rate of T0 transport into the tumor compartment calculated based on T cell transmigration rate (Zhu 1996, PMID: 8706023) and T cell adhesion density (Zhu 1996, PMID: 8706023) [1/day*dm^3]")  # Table S4 (q_T0_T_in = 5.8e-05 1/(centimeter^3*minute))
    q_nT0_P_in <- fixed(254.016); label("rate of nT0 transport into the peripheral compartment calculated based on T cell transmigration rate (Zhu 1996, PMID: 8706023) and T cell adhesion density (Zhu 1996, PMID: 8706023) [1/day]")  # Table S4 (q_nT0_P_in = 0.1764 1/minute)
    q_nT0_P_out <- fixed(5.1); label("rate of nT0 transport out of the peripheral compartment calculated based on steady-state naive T cell density in healthy individuals (Autissier 2010, PMID: 20099249) [1/day]")  # Table S4 (q_nT0_P_out = 5.1 1/day)
    Q_nT0_thym <- fixed(70000000); label("Thymic output of nT0 (Bains 2009, PMID: 19179300; Yeh 2017, PMID: 28912973) [cell/day]")  # Table S4 (Q_nT0_thym = 7e+07 cell/day)
    k_nT0_pro <- fixed(320000000); label("rate of nT0 proliferation (Braber 2012, PMID: 22365666) [cell/day]")  # Table S4 (k_nT0_pro = 3.2e+08 cell/day)
    K_nT0_pro <- fixed(1000000000); label("half-maximal peripheral proliferation of of nT0 (Braber 2012, PMID: 22365666) [cell]")  # Table S4 (K_nT0_pro = 1e+09 cell)
    k_nT0_death <- fixed(0.002); label("rate of nT0 death (Braber 2012, PMID: 22365666) [1/day]")  # Table S4 (k_nT0_death = 0.002 1/day)
    k_IL2_deg <- fixed(288); label("rate of IL2 degradation (Lotze 1985, PMID: 3871099) [1/day]")  # Table S4 (k_IL2_deg = 0.2 1/minute)
    k_IL2_cons <- fixed(0.000144); label("rate of IL2 consumption by T cells (Lotze 1985, PMID: 3871099) [nmol/day*cell]")  # Table S4 (k_IL2_cons = 6e-06 nanomole/cell/h)
    k_IL2_sec <- fixed(0.00072); label("rate of IL2 secretion from T cells (Han 2012, PMID: 22160692; Thurley 2015, PMID: 25923703; Liu 2001) [nmol/day*cell]")  # Table S4 (k_IL2_sec = 3e-05 nanomole/cell/h)
    IL2_50 <- fixed(0.32); label("T cell activation half-maximal IL2 concentration (Marchingo 2014, PMID: 25430770) [nmol/dm^3]")  # Table S4 (IL2_50 = 0.32 nanomolarity)
    IL2_50_Treg <- fixed(0.32); label("Treg activation half-maximal IL2 concentration (Wang and Smith 1987, PMID: 3116143,doi: 10.1084/jem.166.4.1055) [nmol/dm^3]")  # Table S4 (IL2_50_Treg = 0.32 nanomolarity)
    N0 <- fixed(2); label("numer of activated T cell generation by TCR signaling only(doi: 10.1126/science.1260044) [1]")  # Table S4 (N0 = 2 dimensionless)
    N_costim <- fixed(3); label("numer of activated T cell generation by co-stimulatory signaling only (Marchingo 2014, PMID: 25430770) [1]")  # Table S4 (N_costim = 3 dimensionless)
    N_IL2_CD8 <- fixed(11); label("maximum number of activated T cell generations due to IL2 (Marchingo 2014, PMID: 25430770) [1]")  # Table S4 (N_IL2_CD8 = 11 dimensionless)
    N_IL2_CD4 <- fixed(8); label("maximum number of activated T cell generations due to IL2 (Marchingo 2014, PMID: 25430770) [1]")  # Table S4 (N_IL2_CD4 = 8 dimensionless)
    k_Treg <- fixed(0.1); label("Rate of T cell death by Tregs (estimated) [1/day]")  # Table S4 (k_Treg = 0.1 1/day)
    div_T1 <- fixed(1110000); label("T Cell Diversity (Robins 2009, PMID: 19706884) [1]")  # Table S4 (div_T1 = 1.11e+06 dimensionless)
    n_T1_clones <- fixed(92); label("Number of T cell clones (Jafarnejad 2019 PMID: 31236847)) [1]")  # Table S4 (n_T1_clones = 92 dimensionless)
    q_nT1_LN_in <- fixed(0.076); label("Rate of naive T cell transport into the LN calculated based on the steady-state naive T cell level in healthy individuals (Autissier 2010, PMID: 20099249) [1/day]")  # Table S4 (q_nT1_LN_in = 0.076 1/day)
    q_T1_LN_out <- fixed(24); label("Rate of activated T cell transport out of the LN (De Boer and Perelson 1995, PMID: 7475092) [1/day]")  # Table S4 (q_T1_LN_out = 24 1/day)
    q_nT1_LN_out <- fixed(1.8); label("Rate of naive T cell transport out of the LN (Mandl 2012, PMID: 23071319) [1/day]")  # Table S4 (q_nT1_LN_out = 1.8 1/day)
    k_T1_act <- fixed(23); label("T1 activation rate (De Boer and Perelson 1995, PMID: 7475092; Robertson-Tessi 2012, PMID: 22051568) [1/day]")  # Table S4 (k_T1_act = 23 1/day)
    k_T1_pro <- fixed(1); label("T1 proliferation rate (Marchingo 2014, PMID: 25430770) [1/day]")  # Table S4 (k_T1_pro = 1 1/day)
    k_T1_death <- fixed(0.01); label("T1 death rate (De Boer and Perelson 1995, PMID: 7475092) [1/day]")  # Table S4 (k_T1_death = 0.01 1/day)
    q_T1_P_in <- fixed(3507.84); label("rate of T1 transport into the peripheral compartment calculated based on T cell transmigration rate (Zhu 1996, PMID: 8706023) and T cell adhesion density (Zhu 1996, PMID: 8706023) [1/day]")  # Table S4 (q_T1_P_in = 2.436 1/minute)
    q_T1_P_out <- fixed(24); label("rate of T1 transport out of the peripheral compartment (De Boer and Perelson 1995, PMID: 7475092) [1/day]")  # Table S4 (q_T1_P_out = 24 1/day)
    q_T1_T_in <- fixed(83.52); label("rate of T1 transport into the tumor compartment calculated based on T cell transmigration rate (Zhu 1996, PMID: 8706023) and T cell adhesion density (Zhu 1996, PMID: 8706023) [1/day*dm^3]")  # Table S4 (q_T1_T_in = 5.8e-05 1/(centimeter^3*minute))
    q_nT1_P_in <- fixed(254.016); label("rate of nT1 transport into the peripheral compartment calculated based on T cell transmigration rate (Zhu 1996, PMID: 8706023) and T cell adhesion density (Zhu 1996, PMID: 8706023) [1/day]")  # Table S4 (q_nT1_P_in = 0.1764 1/minute)
    q_nT1_P_out <- fixed(5.1); label("rate of nT1 transport out of the peripheral compartment calculated based on steady-state naive T cell density in healthy individuals (Autissier 2010, PMID: 20099249) [1/day]")  # Table S4 (q_nT1_P_out = 5.1 1/day)
    Q_nT1_thym <- fixed(35000000); label("Thymic output of nT1 (Bains 2009, PMID: 19179300; Ye 2002, PMID: 11994448) [cell/day]")  # Table S4 (Q_nT1_thym = 3.5e+07 cell/day)
    k_nT1_pro <- fixed(320000000); label("rate of nT1 proliferation (Braber 2012, PMID: 22365666) [cell/day]")  # Table S4 (k_nT1_pro = 3.2e+08 cell/day)
    K_nT1_pro <- fixed(1000000000); label("half-maximal peripheral proliferation of of nT1 (Braber 2012, PMID: 22365666) [cell]")  # Table S4 (K_nT1_pro = 1e+09 cell)
    k_nT1_death <- fixed(0.002); label("rate of nT1 death (Braber 2012, PMID: 22365666) [1/day]")  # Table S4 (k_nT1_death = 0.002 1/day)
    k_T1 <- fixed(0.1); label("Rate of T cell exhaustion by cancer cells (estimated) [1/day]")  # Table S4 (k_T1 = 0.1 1/day)
    k_C_T1 <- fixed(4); label("Rate of cancer cell death by T cells (Robertson-Tessi 2012, PMID: 22051568) [1/day]")  # Table S4 (k_C_T1 = 4 1/day)
    k_APC_mat <- fixed(1.5); label("Maximum rate of APC maturation (Chen 2014, PMID: 25184733) [1/day]")  # Table S4 (k_APC_mat = 1.5 1/day)
    k_APC_mig <- fixed(4); label("Rate of APC migration (Russo Halin 2016, PMID: 26876174) [1/day]")  # Table S4 (k_APC_mig = 4 1/day)
    k_APC_death <- fixed(0.01); label("Rate of APC death (Marino and Kirschner 2004, PMID: 15038983) [1/day]")  # Table S4 (k_APC_death = 0.01 1/day)
    k_mAPC_death <- fixed(0.02); label("Rate of mAPC death(Marino and Kirschner 2004, PMID: 15038983) [1/day]")  # Table S4 (k_mAPC_death = 0.02 1/day)
    APC0_T <- fixed(400000000); label("Steady state density of CD103+ APCs in the tumor, calculated to fit basal APC equal to 26000/cm^3 from (Lavin 2017, PMID: 28475900) [cell/dm^3]")  # Table S4 (APC0_T = 400000 cell/milliliter)
    APC0_LN <- fixed(1200000000); label("Steady state density of all APC subtypes in the LN (Catron 2004, PMID: 15357945)- divided by 10 to only take CD103+ based on (Lavin 2017, PMID: 28475900) [cell/dm^3]")  # Table S4 (APC0_LN = 1.2e+06 cell/milliliter)
    k_c <- fixed(2); label("Cytokine rate constant (Chen 2014, PMID: 25184733) [1/day]")  # Table S4 (k_c = 2 1/day)
    c0 <- fixed(1); label("Baseline cytokine concentration (Chen 2014, PMID: 25184733) [nmol/dm^3]")  # Table S4 (c0 = 1e-09 molarity)
    c50 <- fixed(1); label("Cytokine concentration for half-maximal APC maturation (Chen 2014, PMID: 25184733) [nmol/dm^3]")  # Table S4 (c50 = 1e-09 molarity)
    DAMPs <- fixed(1.34e-05); label("Concentration of cytokines released per dying cancer cell (Milo 2013, PMID: 24114984, Ponomarenko 2016, PMID: 27298622) [nmol/cell]")  # Table S4 (DAMPs = 1.34e-14 mole/cell)
    n_sites_APC <- fixed(10); label("Maxium number of T Cells an APC can interact with (De Boer and Perelson 1995, PMID: 7475092) [1]")  # Table S4 (n_sites_APC = 10 dimensionless)
    kin <- fixed(14.4); label("Rate of MHC internalization (Chen 2014, PMID: 25184733) [1/day]")  # Table S4 (kin = 14.4 1/day)
    kout <- fixed(28.8); label("Rate of MHC externalization (Chen 2014, PMID: 25184733) [1/day]")  # Table S4 (kout = 28.8 1/day)
    k_P0_up <- fixed(14.4); label("Rate of antigen uptake by APCs (Chen 2014, PMID: 25184733) [1/day*cell]")  # Table S4 (k_P0_up = 14.4 1/day/cell)
    k_xP0_deg <- fixed(2); label("Rate of extracellular antigen degradation (Palsson 2013, PMID: 24074340) [1/day]")  # Table S4 (k_xP0_deg = 2 1/day)
    k_P0_deg <- fixed(17.28); label("Rate of endosomal antigen degradation (Chen 2014, PMID: 25184733) [1/day]")  # Table S4 (k_P0_deg = 17.28 1/day)
    k_p0_deg <- fixed(144); label("Rate of endosomal epitope degradation (Chen 2014, PMID: 25184733) [1/day]")  # Table S4 (k_p0_deg = 144 1/day)
    k_P0_on <- fixed(0.000144); label("Rate of antigen-MHC binding (Agrawal and Linderman 1996, PMID: 8944895) [dm^3/day*nmol]")  # Table S4 (k_P0_on = 144000 1/day/molarity)
    k_P0_d1 <- fixed(100); label("Antigen-MHC kd (Stone 2015, PMID: 25618219) [nmol/dm^3]")  # Table S4 (k_P0_d1 = 1e-07 molarity)
    p0_50 <- fixed(4.39296049517e-10); label("calculated based on the number of molecules for half-maximal T cell activation (Kimachi 1997, PMID: 9464819) and synapse surface area [nmol/dm^2]")  # Table S4 (p0_50 = 2.6455e-05 molecule/micrometer^2)
    P0_C1 <- fixed(0.00054); label("Concentration of P0 in C1 (estimated by 2.7e6 proteins/um^3 in a cancer cell (Milo 2013, PMID: 24114984) * cancer cell volume / Avogadro constant) [nmol/cell]")  # Table S4 (P0_C1 = 5.4e-13 mole/cell)
    P0_C2 <- fixed(0.00054); label("Concentration of P0 in C2 (estimated by 2.7e6 proteins/um^3 in a cancer cell (Milo 2013, PMID: 24114984) * cancer cell volume / Avogadro constant) [nmol/cell]")  # Table S4 (P0_C2 = 5.4e-13 mole/cell)
    A_Tcell <- fixed(1.5131041193e-08); label("Surface area of the T cell calculated based on the average T cell diameter (Chapman 1981, PMID: 6975780) [dm^2]")  # Table S4 (A_Tcell = 151.31 micrometer^2)
    A_cell <- fixed(8.97270277792e-08); label("Surface area of the Cancer cell calculated based on the average cancer cell diameter (Jafarnejad 2019 PMID: 31236847) [dm^2]")  # Table S4 (A_cell = 897.27 micrometer^2)
    A_APC <- fixed(9e-08); label("Surface area of the APC (Agrawal and Linderman 1996, PMID: 8944895) [dm^2]")  # Table S4 (A_APC = 900 micrometer^2)
    k_M1p0_TCR_on <- fixed(5203129616.64); label("Rate of TCR binding to MHC-peptide complex (Lever 2014, PMID: 25145757) [dm^2/day*nmol]")  # Table S4 (k_M1p0_TCR_on = 1 1/(second*molecule/micrometer^2))
    k_M1p0_TCR_off <- fixed(86400); label("Rate of TCR unbinding from MHC-peptide complex (Lever 2014, PMID: 25145757) [1/day]")  # Table S4 (k_M1p0_TCR_off = 1 1/second)
    k_M1p0_TCR_p <- fixed(86400); label("Rate of MHC-peptide-TCR complex modification (Lever 2014, PMID: 25145757) [1/day]")  # Table S4 (k_M1p0_TCR_p = 1 1/second)
    phi_M1p0_TCR <- fixed(7776); label("Rate of MHC-peptide-TCR complex with maximal modification that leads to non-signaling(Lever 2014, PMID: 25145757) [1/day]")  # Table S4 (phi_M1p0_TCR = 0.09 1/second)
    N_M1p0_TCR <- fixed(10); label("Number of modifcation steps for MHC-peptide-TCR complex (Lever 2014, PMID: 25145757) [1]")  # Table S4 (N_M1p0_TCR = 10 dimensionless)
    TCR_p0_tot <- fixed(0.00172385676137); label("Total number of TCR molecules per naive T cell calculated based on the number of molecules for TCR (Lever 2014, PMID: 25145757) and T cell surface area [nmol/dm^2]")  # Table S4 (TCR_p0_tot = 103.813 molecule/micrometer^2)
    k_P1_up <- fixed(14.4); label("Rate of antigen uptake by APCs (Chen 2014, PMID: 25184733) [1/day*cell]")  # Table S4 (k_P1_up = 14.4 1/day/cell)
    k_xP1_deg <- fixed(2); label("Rate of extracellular antigen degradation (Palsson 2013, PMID: 24074340) [1/day]")  # Table S4 (k_xP1_deg = 2 1/day)
    k_P1_deg <- fixed(17.28); label("Rate of endosomal antigen degradation (Chen 2014, PMID: 25184733) [1/day]")  # Table S4 (k_P1_deg = 17.28 1/day)
    k_p1_deg <- fixed(144); label("Rate of endosomal epitope degradation (Chen 2014, PMID: 25184733) [1/day]")  # Table S4 (k_p1_deg = 144 1/day)
    k_P1_on <- fixed(0.000144); label("Rate of antigen-MHC binding (Agrawal and Linderman 1996, PMID: 8944895) [dm^3/day*nmol]")  # Table S4 (k_P1_on = 144000 1/day/molarity)
    k_P1_d1 <- fixed(10); label("Antigen-MHC kd (Stone 2015, PMID: 25618219) [nmol/dm^3]")  # Table S4 (k_P1_d1 = 1e-08 molarity)
    p1_50 <- fixed(4.39296049517e-10); label("calculated based on the number of molecules for half-maximal T cell activation (Kimachi 1997, PMID: 9464819) and synapse surface area [nmol/dm^2]")  # Table S4 (p1_50 = 2.6455e-05 molecule/micrometer^2)
    P1_C1 <- fixed(0.00054); label("Concentration of P1 in C1 (estimated by 2.7e6 proteins/um^3 in a cancer cell (Milo 2013, PMID: 24114984) * cancer cell volume / Avogadro constant) [nmol/cell]")  # Table S4 (P1_C1 = 5.4e-13 mole/cell)
    P1_C2 <- fixed(0.00054); label("Concentration of P1 in C2 (estimated by 2.7e6 proteins/um^3 in a cancer cell (Milo 2013, PMID: 24114984) * cancer cell volume / Avogadro constant) [nmol/cell]")  # Table S4 (P1_C2 = 5.4e-13 mole/cell)
    k_M1p1_TCR_on <- fixed(5203129616.64); label("Rate of TCR binding to MHC-peptide complex (Lever 2014, PMID: 25145757) [dm^2/day*nmol]")  # Table S4 (k_M1p1_TCR_on = 1 1/(second*molecule/micrometer^2))
    k_M1p1_TCR_off <- fixed(86400); label("Rate of TCR unbinding from MHC-peptide complex (Lever 2014, PMID: 25145757) [1/day]")  # Table S4 (k_M1p1_TCR_off = 1 1/second)
    k_M1p1_TCR_p <- fixed(86400); label("Rate of MHC-peptide-TCR complex modification (Lever 2014, PMID: 25145757) [1/day]")  # Table S4 (k_M1p1_TCR_p = 1 1/second)
    phi_M1p1_TCR <- fixed(7776); label("Rate of MHC-peptide-TCR complex with maximal modification that leads to non-signaling(Lever 2014, PMID: 25145757) [1/day]")  # Table S4 (phi_M1p1_TCR = 0.09 1/second)
    N_M1p1_TCR <- fixed(10); label("Number of modifcation steps for MHC-peptide-TCR complex (Lever 2014, PMID: 25145757) [1]")  # Table S4 (N_M1p1_TCR = 10 dimensionless)
    TCR_p1_tot <- fixed(0.00172385676137); label("Total number of TCR molecules per naive T cell calculated based on the number of molecules for TCR (Lever 2014, PMID: 25145757) and T cell surface area [nmol/dm^2]")  # Table S4 (TCR_p1_tot = 103.813 molecule/micrometer^2)
    kon_PD1_PDL1 <- fixed(504000000); label("kon of PD1-PDL1 binding kon for PD1-PDL1 (Cheng 2013, PMID: 23417675) [dm^2/day*nmol]")  # Table S4 (kon_PD1_PDL1 = 0.0583333 1/(micromolarity*nanometer*second))
    q_P_nivo <- fixed(0.76032); label("Volumetric flow rate of nivo between central and peripheral compartment [dm^3/day]")  # Table S4 (q_P_nivo = 8.8e-06 liter/second)
    q_T_nivo <- fixed(0.07475328); label("Volumetric flow rate of nivo between central and tumor compartment [dm^3/day]")  # Table S4 (q_T_nivo = 0.0008652 milliliter/second)
    q_LN_nivo <- fixed(0.0002808); label("Volumetric flow rate of nivo between central and TDLN compartment [dm^3/day]")  # Table S4 (q_LN_nivo = 3.25e-06 milliliter/second)
    q_LD_nivo <- fixed(2.16); label("Rate of lymphatic drainage of nivo from TDLN to central compartment [1/day]")  # Table S4 (q_LD_nivo = 0.0015 1/minute)
    k_cl_nivo <- fixed(0.37); label("Clearance rate of nivo from central compartment [dm^3/day]")  # Table S4 (k_cl_nivo = 0.37 liter/day)
    gamma_C_nivo <- fixed(0.61); label("Volume fraction of interstitial space available to nivo in central compartment [1]")  # Table S4 (gamma_C_nivo = 0.61 dimensionless)
    gamma_P_nivo <- fixed(0.0452); label("Volume fraction of interstitial space available to nivo in peripheral compartment [1]")  # Table S4 (gamma_P_nivo = 0.0452 dimensionless)
    gamma_T_nivo <- fixed(0.522); label("Volume fraction of interstitial space available to nivo in tumor compartment [1]")  # Table S4 (gamma_T_nivo = 0.522 dimensionless)
    gamma_LN_nivo <- fixed(0.2); label("Volume fraction of interstitial space available to nivo in TDLN compartment [1]")  # Table S4 (gamma_LN_nivo = 0.2 dimensionless)
    q_P_atezo <- fixed(0.75168); label("Volumetric flow rate of atezo between central and peripheral compartment [dm^3/day]")  # Table S4 (q_P_atezo = 8.7e-06 liter/second)
    q_T_atezo <- fixed(0.00736128); label("Volumetric flow rate of atezo between central and tumor compartment [dm^3/day]")  # Table S4 (q_T_atezo = 8.52e-05 milliliter/second)
    q_LN_atezo <- fixed(0.0002808); label("Volumetric flow rate of atezo between central and TDLN compartment [dm^3/day]")  # Table S4 (q_LN_atezo = 3.25e-06 milliliter/second)
    q_LD_atezo <- fixed(2.16); label("Rate of lymphatic drainage of atezo from TDLN to central compartment [1/day]")  # Table S4 (q_LD_atezo = 0.0015 1/minute)
    k_cl_atezo <- fixed(0.324); label("Clearance rate of atezo from central compartment [dm^3/day]")  # Table S4 (k_cl_atezo = 0.324 liter/day)
    gamma_C_atezo <- fixed(0.61); label("Volume fraction of interstitial space available to atezo in central compartment [1]")  # Table S4 (gamma_C_atezo = 0.61 dimensionless)
    gamma_P_atezo <- fixed(0.068); label("Volume fraction of interstitial space available to atezo in peripheral compartment [1]")  # Table S4 (gamma_P_atezo = 0.068 dimensionless)
    gamma_T_atezo <- fixed(0.522); label("Volume fraction of interstitial space available to atezo in tumor compartment [1]")  # Table S4 (gamma_T_atezo = 0.522 dimensionless)
    gamma_LN_atezo <- fixed(0.2); label("Volume fraction of interstitial space available to atezo in TDLN compartment [1]")  # Table S4 (gamma_LN_atezo = 0.2 dimensionless)
    q_P_ipi <- fixed(0.9936); label("Volumetric flow rate of ipi between central and peripheral compartment [dm^3/day]")  # Table S4 (q_P_ipi = 1.15e-05 liter/second)
    q_T_ipi <- fixed(0.00736128); label("Volumetric flow rate of ipi between central and tumor compartment [dm^3/day]")  # Table S4 (q_T_ipi = 8.52e-05 milliliter/second)
    q_LN_ipi <- fixed(0.0002808); label("Volumetric flow rate of ipi between central and TDLN compartment [dm^3/day]")  # Table S4 (q_LN_ipi = 3.25e-06 milliliter/second)
    q_LD_ipi <- fixed(2.16); label("Rate of lymphatic drainage of ipi from TDLN to central compartment [1/day]")  # Table S4 (q_LD_ipi = 0.0015 1/minute)
    k_cl_ipi <- fixed(0.524); label("Clearance rate of ipi from central compartment [dm^3/day]")  # Table S4 (k_cl_ipi = 0.524 liter/day)
    gamma_C_ipi <- fixed(0.686); label("Volume fraction of interstitial space available to ipi in central compartment [1]")  # Table S4 (gamma_C_ipi = 0.686 dimensionless)
    gamma_P_ipi <- fixed(0.0496); label("Volume fraction of interstitial space available to ipi in peripheral compartment [1]")  # Table S4 (gamma_P_ipi = 0.0496 dimensionless)
    gamma_T_ipi <- fixed(0.522); label("Volume fraction of interstitial space available to ipi in tumor compartment [1]")  # Table S4 (gamma_T_ipi = 0.522 dimensionless)
    gamma_LN_ipi <- fixed(0.2); label("Volume fraction of interstitial space available to ipi in TDLN compartment [1]")  # Table S4 (gamma_LN_ipi = 0.2 dimensionless)
    k_out_PDL1 <- fixed(8.30269533587e-11); label("Expression rate of PDL1 on tumor cells (Mimura 2018, PMID: 29034543) [nmol/day]")  # Table S4 (k_out_PDL1 = 50000 molecule/day)
    k_in_PDL1 <- fixed(1); label("Degradation rate of PDL1 on tumor cells (Hsu 2018, PMID: 30442814) [1/day]")  # Table S4 (k_in_PDL1 = 1 1/day)
    r_PDL1_IFNg <- fixed(6); label("Number of folds increase of PDL1 expression by IFNg (Shin 2017, PMID: 27903500) [1]")  # Table S4 (r_PDL1_IFNg = 6 dimensionless)
    kon_PD1_PDL2 <- fixed(662400000); label("kon of PD1-PDL2 binding (Cheng 2013, PMID: 23417675) [dm^2/day*nmol]")  # Table S4 (kon_PD1_PDL2 = 0.0766667 1/(micromolarity*nanometer*second))
    kon_PD1_nivo <- fixed(112.32); label("kon of PD1-nivolumab binding (Wang 2014, PMID: 24872026) [dm^3/day*nmol]")  # Table S4 (kon_PD1_nivo = 1.3e+06 1/(molarity*second))
    kon_PDL1_atezo <- fixed(37.152); label("kon of PDL1-atezolizumab binding (Wang 2014, PMID: 24872026) [dm^3/day*nmol]")  # Table S4 (kon_PDL1_atezo = 430000 1/(molarity*second))
    kon_CD28_CD80 <- fixed(1152000000); label("kon of CD28-CD80 binding (Jansson 2005, PMID: 16034096) [dm^2/day*nmol]")  # Table S4 (kon_CD28_CD80 = 0.133333 1/(micromolarity*nanometer*second))
    kon_CD28_CD86 <- fixed(4032000000); label("kon of CD28-CD86 binding (Jansson 2005, PMID: 16034096) [dm^2/day*nmol]")  # Table S4 (kon_CD28_CD86 = 0.466667 1/(micromolarity*nanometer*second))
    kon_CTLA4_CD80 <- fixed(6336000000); label("kon of CTLA4-CD80 binding (Jansson 2005, PMID: 16034096) [dm^2/day*nmol]")  # Table S4 (kon_CTLA4_CD80 = 0.733333 1/(micromolarity*nanometer*second))
    kon_CTLA4_CD86 <- fixed(5760000000); label("kon of CTLA4-CD86 binding (Jansson 2005, PMID: 16034096) [dm^2/day*nmol]")  # Table S4 (kon_CTLA4_CD86 = 0.666667 1/(micromolarity*nanometer*second))
    kon_CD80_PDL1 <- fixed(9100800000); label("kon of CD80-PDL1 binding (Jansson 2005, PMID: 16034096) [dm^2/day*nmol]")  # Table S4 (kon_CD80_PDL1 = 1.05333 1/(micromolarity*nanometer*second))
    kon_CTLA4_ipi <- fixed(33.0912); label("kon of CTLA4-ipilimumab binding (Wang 2014, PMID: 24872026) [dm^3/day*nmol]")  # Table S4 (kon_CTLA4_ipi = 383000 1/(molarity*second))
    koff_PD1_PDL1 <- fixed(123984); label("koff of PD1-PDL1 binding calculated based on the measured kd and kon (Cheng 2013, PMID: 23417675) [1/day]")  # Table S4 (koff_PD1_PDL1 = 1.435 1/second)
    koff_PD1_PDL2 <- fixed(45705.6); label("koff of PD1-PDL2 binding calculated based on the measured kd and kon (Cheng 2013, PMID: 23417675) [1/day]")  # Table S4 (koff_PD1_PDL2 = 0.529 1/second)
    koff_PD1_nivo <- fixed(292.032); label("koff of PD1-nivolumab binding calculated based on the measured kd and kon (Wang 2014, PMID: 24872026) [1/day]")  # Table S4 (koff_PD1_nivo = 0.00338 1/second)
    koff_PDL1_atezo <- fixed(14.8608); label("koff of PDL1-atezolizumab binding calculated based on the measured kd and kon (Wang 2014, PMID: 24872026) [1/day]")  # Table S4 (koff_PDL1_atezo = 0.000172 1/second)
    koff_CD28_CD80 <- fixed(138240); label("koff of CD28-CD80 binding calculated based on the measured kd and kon (Jansson 2005, PMID: 16034096) [1/day]")  # Table S4 (koff_CD28_CD80 = 1.6 1/second)
    koff_CD28_CD86 <- fixed(2419200); label("koff of CD28-CD86 binding calculated based on the measured kd and kon (Jansson 2005, PMID: 16034096) [1/day]")  # Table S4 (koff_CD28_CD86 = 28 1/second)
    koff_CTLA4_CD80 <- fixed(38016); label("koff of CTLA4-CD80 binding calculated based on the measured kd and kon (Jansson 2005, PMID: 16034096) [1/day]")  # Table S4 (koff_CTLA4_CD80 = 0.44 1/second)
    koff_CTLA4_CD86 <- fixed(449280); label("koff of CTLA4-CD86 binding calculated based on the measured kd and kon (Jansson 2005, PMID: 16034096) [1/day]")  # Table S4 (koff_CTLA4_CD86 = 5.2 1/second)
    koff_CD80_PDL1 <- fixed(518745.6); label("koff of CD80-PDL1 binding calculated based on the measured kd and kon (Jansson 2005, PMID: 16034096) [1/day]")  # Table S4 (koff_CD80_PDL1 = 6.004 1/second)
    koff_CTLA4_ipi <- fixed(602.25984); label("koff of CTLA4-ipilimumab binding calculated based on the measured kd and kon (Wang 2014, PMID: 24872026) [1/day]")  # Table S4 (koff_CTLA4_ipi = 0.0069706 1/second)
    Chi_PD1_nivo <- fixed(333333330.953); label("Antibody cross-arm binding efficiency that also includes the conversion of kon from 3D to 2D (estimated based on Wang 2014, PMID: 24872026) [1/dm]")  # Table S4 (Chi_PD1_nivo = 3.33333 1/nanometer)
    Chi_PDL1_atezo <- fixed(3333333309.53); label("Antibody cross-arm binding efficiency that also includes the conversion of kon from 3D to 2D (Wang 2014, PMID: 24872026) [1/dm]")  # Table S4 (Chi_PDL1_atezo = 33.3333 1/nanometer)
    Chi_CTLA4_ipi <- fixed(333333330.953); label("Antibody cross-arm binding efficiency that also includes the conversion of kon from 3D to 2D (estimated) [1/dm]")  # Table S4 (Chi_CTLA4_ipi = 3.33333 1/nanometer)
    PD1_50 <- fixed(0.000830269533587); label("PD1/PDL1 concentration for half-maximal T cell inactivation (estimated) [nmol/dm^2]")  # Table S4 (PD1_50 = 50 molecule/micrometer^2)
    n_PD1 <- fixed(2); label("Hill coefficient for PD1/PDL1 half-maximal T cell inactivation (estimated) [1]")  # Table S4 (n_PD1 = 2 dimensionless)
    CD28_CD8X_50 <- fixed(0.00332107813435); label("CD28-CD80/CD28-CD86 concentration for half-maximal T cell co-estimulation (estimated) [nmol/dm^2]")  # Table S4 (CD28_CD8X_50 = 200 molecule/micrometer^2)
    n_CD28_CD8X <- fixed(2); label("Hill coefficient for CD28-CD80/CD28-CD86 half-maximal T cell co-estimulation (estimated) [1]")  # Table S4 (n_CD28_CD8X = 2 dimensionless)
    T_PD1_total <- fixed(1.02953422165e-10); label("concentration of PD1 on T cells (Cheng 2014, PMID: 23417675; Mkrtichyan 2012, PMID: 22837483) [nmol]")  # Table S4 (T_PD1_total = 62000 molecule)
    T_CD28_total <- fixed(3.0553918836e-10); label("concentration of CD28 on T cells (Jansson 2005, PMID: 16034096) [nmol]")  # Table S4 (T_CD28_total = 184000 molecule)
    T_CTLA4_syn <- fixed(1.32843125374e-11); label("concentration of CTLA4 on T cells (Jansson 2005, PMID: 16034096) [nmol]")  # Table S4 (T_CTLA4_syn = 8000 molecule)
    T_PDL1_total <- fixed(3.08860266494e-10); label("concentration of PDL1 on T cells (Cheng 2014, PMID: 23417675; Mkrtichyan 2012, PMID: 22837483) [nmol]")  # Table S4 (T_PDL1_total = 186000 molecule)
    C1_PDL1_base <- fixed(2.98897032091e-10); label("baseline number of PDL1 molecules per C1 cell (Cheng 2014, PMID: 23417675; Mkrtichyan 2012, PMID: 22837483; Gatalica 2014, PMID: 25392179) [nmol]")  # Table S4 (C1_PDL1_base = 180000 molecule)
    r_PDL2C1 <- fixed(0.07); label("PDL2/PDL1 molecule ratio on C1cell (Cheng 2014, PMID: 23417675) [1]")  # Table S4 (r_PDL2C1 = 0.07 dimensionless)
    C1_CD80_total <- fixed(6.6421562687e-11); label("number of CD80 molecules per C1 cell (Jansson 2005, PMID: 16034096) [nmol]")  # Table S4 (C1_CD80_total = 40000 molecule)
    C1_CD86_total <- fixed(1.42806359777e-09); label("number of CD86 molecules per C1 cell (Jansson 2005, PMID: 16034096) [nmol]")  # Table S4 (C1_CD86_total = 860000 molecule)
    APC_PDL1_base <- fixed(4.42810417913e-10); label("baseline number of PDL1 molecules per APC cell (Cheng 2014, PMID: 23417675; Mkrtichyan 2012, PMID: 22837483; Gatalica 2014, PMID: 25392179) [nmol]")  # Table S4 (APC_PDL1_base = 266667 molecule)
    r_PDL2APC <- fixed(0.07); label("PDL2/PDL1 molecule ratio on APCcell (Cheng 2014, PMID: 23417675) [1]")  # Table S4 (r_PDL2APC = 0.07 dimensionless)
    APC_CD80_total <- fixed(6.6421562687e-11); label("number of CD80 molecules per APC cell (Jansson 2005, PMID: 16034096) [nmol]")  # Table S4 (APC_CD80_total = 40000 molecule)
    APC_CD86_total <- fixed(1.42806359777e-09); label("number of CD86 molecules per APC cell (Jansson 2005, PMID: 16034096) [nmol]")  # Table S4 (APC_CD86_total = 860000 molecule)
    Treg_CTLA4_tot <- fixed(8.30269533587e-12); label("Total number of CTLA4 on Treg cells (Jansson 2005, PMID: 16034096) [nmol]")  # Table S4 (Treg_CTLA4_tot = 5000 molecule)
    Treg_CTLA4_50 <- fixed(1.66053906717e-12); label("CTLA4 occupancy for half-maximal Treg inactivation by macrophages (estimated) [nmol]")  # Table S4 (Treg_CTLA4_50 = 1000 molecule)
    n_Treg_CTLA4 <- fixed(2); label("CTLA4 occupancy Hill coefficient for Treg inactivation by macrophages (estimated) [1]")  # Table S4 (n_Treg_CTLA4 = 2 dimensionless)
    k_CTLA4_ADCC <- fixed(0.1); label("Anti-CTLA4 ADCC (antibody-dependent cellular cytotoxicity) rate of Treg (Richards 2008, PMID: 18723496) [1/day]")  # Table S4 (k_CTLA4_ADCC = 0.1 1/day)
    k_rec_MDSC <- fixed(1.2); label("Rate of MDSC recruitment into the tumor (Lai 2018, PMID: 29735668; Huang 2006, PMID: 17257744) [1/day]")  # Table S4 (k_rec_MDSC = 1.2 1/day)
    kd_MDSC <- fixed(0.015); label("Rate of MDSC death (Lai 2018, PMID: 29735668) [1/day]")  # Table S4 (kd_MDSC = 0.015 1/day)
    IC50_ENT_C <- fixed(374); label("ENT concentration for half-maximal tumor cell death (Bouchain 2002, PMID: 12593661; Lee 2001, PMID: 11221885) [nmol/dm^3]")  # Table S4 (IC50_ENT_C = 3.74e-07 molarity)
    k_deg_CCL2 <- fixed(1.44); label("rate of CCL2 degradation (Tanimoto 2007, PMID: 18089573) [1/day]")  # Table S4 (k_deg_CCL2 = 0.06 1/h)
    k_deg_NO <- fixed(135); label("rate of NO degradation (Hakim 1996, PMID: 8953625) [1/day]")  # Table S4 (k_deg_NO = 135 1/day)
    k_deg_ArgI <- fixed(0.173); label("rate of ArgI degradation (Schimke 1964, PMID: 14257612) [1/day]")  # Table S4 (k_deg_ArgI = 0.173 1/day)
    k_sec_CCL2 <- fixed(1.42e-10); label("rate of CCL2 secretion (Huang 2007, PMID: 17257744; Dutta, PMID: 29594759) [nmol/day*cell]")  # Table S4 (k_sec_CCL2 = 1.42e-10 nanomole/cell/day)
    k_sec_NO <- fixed(4.8e-07); label("rate of NO secretion from MDSCs (Serafini 2008, PMID: 18593947) [nmol/day*cell]")  # Table S4 (k_sec_NO = 4.8e-07 nanomole/cell/day)
    k_sec_ArgI <- fixed(1.4e-08); label("rate of ArgI secretion (Serafini 2008, PMID: 18593947) [dm^3*mU/day*cell]")  # Table S4 (k_sec_ArgI = 0.014 (mU*microliter)/cell/day)
    IC50_ENT_NO <- fixed(0.56); label("half-maximal ENT concentration for NO inhibition (Choo 2010, PMID: 20421217) [nmol/dm^3]")  # Table S4 (IC50_ENT_NO = 5.6e-10 molarity)
    IC50_ArgI_CTL <- fixed(61.7); label("rate of ArgI-induced T cell death (Serafini 2008, PMID: 18593947) [mU]")  # Table S4 (IC50_ArgI_CTL = 61.7 mU)
    IC50_NO_CTL <- fixed(0.75); label("rate of NO-induced T cell death (Serafini 2008, PMID: 18593947) [nmol/dm^3]")  # Table S4 (IC50_NO_CTL = 7.5e-10 molarity)
    EC50_CCL2_rec <- fixed(0.5); label("Half-maximal CCL2 level of MDSC recruitment (Ernst 1994, PMID: 8144933) [nmol/dm^3]")  # Table S4 (EC50_CCL2_rec = 5e-10 molarity)
    EC50_ArgI_Treg <- fixed(22.1); label("Half-maximal ArgI level of Treg expansion (estimated based on Serafini 2008, PMID: 18593947) [mU]")  # Table S4 (EC50_ArgI_Treg = 22.1 mU)
    MDSC_max <- fixed(163700000); label("Maximal MDSC density in the tumour (Diaz-Montero 2009, PMID: 18446337) [cell/dm^3]")  # Table S4 (MDSC_max = 163700 cell/milliliter)
    IC50_ENT_CCL2 <- fixed(1.2); label("half-maximal ENT concentration for CCL2 inhibition (Choo 2013, PMID: 24241152) [nmol/dm^3]")  # Table S4 (IC50_ENT_CCL2 = 1.2e-09 molarity)
    k_brec_MDSC <- fixed(0.0021); label("Baseline rate of MDSC migration into the tumor(Huang 2006, PMID: 17257744) [1/day]")  # Table S4 (k_brec_MDSC = 0.0021 1/day)
    IC50_ENT_ArgI <- fixed(500); label("half-maximal ENT concentration for Arg I inhibition (Orillian 2017, PMID: 28698201) [nmol/dm^3]")  # Table S4 (IC50_ENT_ArgI = 5e-07 molarity)
    k_a1_ENT <- fixed(45.6); label("Rate of buccal absorption of ENT into central compartment [1/day]")  # Table S4 (k_a1_ENT = 1.9 1/h)
    k_a2_ENT <- fixed(60); label("Rate of GI absorption of ENT into central compartment [1/day]")  # Table S4 (k_a2_ENT = 2.5 1/h)
    k_cln_ENT <- fixed(2932983.71317); label("Non-linear clearance rate of ENT from central compartment [nmol/day*dm^3]")  # Table S4 (k_cln_ENT = 0.000122208 mole/(liter*hour))
    Kc_ENT <- fixed(9378.10915535); label("Concentration of ENT in central compartment at 50% k_cln [nmol/dm^3]")  # Table S4 (Kc_ENT = 9.37811e-06 mole/liter)
    k_dose2 <- fixed(0.3312); label("Rate of GI absorption of ENT into buccal compartment [1/day]")  # Table S4 (k_dose2 = 0.0138 1/h)
    q_P_ENT <- fixed(27648); label("Volumetric flow rate of ENT between central and peripheral compartment [dm^3/day]")  # Table S4 (q_P_ENT = 320 milliliter/second)
    q_T_ENT <- fixed(276480); label("Volumetric flow rate of ENT between central and tumor compartment [dm^3/day]")  # Table S4 (q_T_ENT = 3200 milliliter/second)
    q_LN_ENT <- fixed(27648); label("Volumetric flow rate of ENT between central and TDLN compartment [dm^3/day]")  # Table S4 (q_LN_ENT = 320 milliliter/second)
    q_LD_ENT <- fixed(2.16); label("Rate of lymphatic drainage of ENT from TDLN to central compartment [1/day]")  # Table S4 (q_LD_ENT = 0.0015 1/minute)
    k_cl_ENT <- fixed(3.648); label("Clearance rate of ENT from central compartment [1/day]")  # Table S4 (k_cl_ENT = 0.152 1/h)
    gamma_C_ENT <- fixed(0.55); label("Volume fraction of interstitial space available to ENT in central compartment [1]")  # Table S4 (gamma_C_ENT = 0.55 dimensionless)
    gamma_P_ENT <- fixed(0.062); label("Volume fraction of interstitial space available to ENT in peripheral compartment [1]")  # Table S4 (gamma_P_ENT = 0.062 dimensionless)
    gamma_T_ENT <- fixed(0.611); label("Volume fraction of interstitial space available to ENT in tumor compartment [1]")  # Table S4 (gamma_T_ENT = 0.611 dimensionless)
    gamma_LN_ENT <- fixed(0.2); label("Volume fraction of interstitial space available to ENT in TDLN compartment [1]")  # Table S4 (gamma_LN_ENT = 0.2 dimensionless)
    Vmcl <- fixed(193680); label("Maximum elimination rate of nab-paclitaxel from the central compartment V1 [ug/day]")  # Table S4 (Vmcl = 8070 microgram/h)
    Kcl <- fixed(40.2); label("Nab-paclitaxel concentration in the central compartment V1 at 50% Vmcl [ug/dm^3]")  # Table S4 (Kcl = 40.2 microgram/liter)
    Q2 <- fixed(998.4); label("Intercompartmental clearance between the central compartment V1 and the second peripheral compartment V3 [dm^3/day]")  # Table S4 (Q2 = 41.6 liter/h)
    Vmt <- fixed(7800000); label("Maximum intercompartmental distribution rate between the central compartment V1 and the first peripheral compartment V2 [ug/day]")  # Table S4 (Vmt = 325000 microgram/h)
    Kt <- fixed(4260); label("Nab-paclitaxel concentration in the central compartment V1 at 50% Vmt [ug/dm^3]")  # Table S4 (Kt = 4260 microgram/liter)
    BSA <- fixed(190); label("Body surface area [dm^2]")  # Table S4 (BSA = 1.9 meter^2)
    r_nabp <- fixed(1.5); label("Tumour to plasma concentration ratio of nab-paclitaxel (Yang 2007, PMID: 17828616; Rajeshkumar 2016, PMID: 27441498) [1]")  # Table S4 (r_nabp = 1.5 dimensionless)
    k_C_nabp <- fixed(1.368); label("Cancer cell killing rate by nab-paclitaxel (Choi 2012, PMID: 23023313) [1/day]")  # Table S4 (k_C_nabp = 0.057 1/h)
    IC50_nabp <- fixed(47); label("Half-maximal nab-paclitaxel concentration for cancer cell killing (Yang 2013, PMID: 23180760) [nmol/dm^3]")  # Table S4 (IC50_nabp = 4.7e-08 molarity)
    MW_nabp <- fixed(0.8539); label("Molecular weight of nab-paclitaxel (calculated) [ug/nmol]")  # Table S4 (MW_nabp = 853.9 gram/mole)
    Kc_cyt <- fixed(80000000); label("Half-Maximal cancer cell number for cytotoxic drug diffusion (estimated based on Kyle 2007, PMID: 17473214) [cell]")  # Table S4 (Kc_cyt = 8e+07 cell)
    k_C_resist <- fixed(5e-05); label("Cancer resistance to nab-paclitaxel (estimated) [1/day]")  # Table S4 (k_C_resist = 5e-05 1/day)
    r_resist <- fixed(100); label("Number of folds increase of nab-paclitaxel EC50 in resistant cancer clones (Nemcova-Furstova 2016, PMID: 27664577) [1]")  # Table S4 (r_resist = 100 dimensionless)
    k_K_g <- fixed(4.93); label("Tumour vasculature growth rate (Desai 2006, PMID: 16489089) [1/day]")  # Table S4 (k_K_g = 4.93 1/day)
    k_K_d <- fixed(0.002); label("Tumour vasculature inhibition rate (Desai 2006, PMID: 16489089) [1/day]")  # Table S4 (k_K_d = 0.002 1/day)
    k_c_nabp <- fixed(1.7e-08); label("Secretion rate of angiogenic factors induced by nab-paclitaxel (Volk 2008, PMID: 18516298) [ug/day*cell]")  # Table S4 (k_c_nabp = 0.017 picogram/cell/day)
    IC50_nabp_vas <- fixed(6.2); label("Half-maximal conc. of nab-paclitaxel on angiogenic factor induction (Volk 2008, PMID: 18516298) [nmol/dm^3]")  # Table S4 (IC50_nabp_vas = 6.2e-09 molarity)
    k_c_vas <- fixed(3.4e-09); label("Secretion rate of angiogenic factors by cancer cells (Volk 2008, PMID: 18516298) [ug/day*cell]")  # Table S4 (k_c_vas = 0.0034 picogram/cell/day)
    k_c_deg <- fixed(11); label("Degradation rate of angiogenic factors (Finley 2011, PMID: 22104283) [1/day]")  # Table S4 (k_c_deg = 11 1/day)
    IC50_vas <- fixed(150); label("Half-maximal conc. of angiogenic factor on tumor capacity growth (Desai 2006, PMID: 16489089) [ug/dm^3]")  # Table S4 (IC50_vas = 150000 picogram/milliliter)
    k_endo_nabp <- fixed(4.272e-07); label("Inhibition rate of maximal tumor capacity by nab-paclitaxel (Mollard 2017, PMID: 28416742) [dm^3/day*nmol]")  # Table S4 (k_endo_nabp = 17.8 1/h/molarity)
    k_Th_act <- fixed(10); label("Rate of T helper cell activation (Robertson-Tessi 2012, PMID: 22051568) [1/day]")  # Table S4 (k_Th_act = 10 1/day)
    k_reg <- fixed(0.022); label("Th differentiation rate to Treg (Robertson-Tessi 2012, PMID: 22051568) [1/day]")  # Table S4 (k_reg = 0.022 1/day)
    k_TGF_Tsec <- fixed(1.2e-10); label("TGF secretion rate by Treg (Liyanage 2002, PMID: 12193750) [nmol/day*cell]")  # Table S4 (k_TGF_Tsec = 1.2e-10 nanomole/cell/day)
    k_TGF_deg <- fixed(14.3); label("TGF degradtion rate (Robertson-Tessi 2012, PMID: 22051568) [1/day]")  # Table S4 (k_TGF_deg = 14.3 1/day)
    TGF_50_reg <- fixed(0.14); label("Half-Maximal TGFb level for differentiation of Th to Treg (Robertson-Tessi 2012, PMID: 22051568) [nmol/dm^3]")  # Table S4 (TGF_50_reg = 1.4e-10 molarity)
    TGF_50_ctl <- fixed(0.28); label("Half-Maximal TGFb level for CD8 T cell inhibition (Robertson-Tessi 2012, PMID: 22051568) [nmol/dm^3]")  # Table S4 (TGF_50_ctl = 2.8e-10 molarity)
    Kc_rec <- fixed(20200000); label("Half-Maximal cancer cell number for T cell recruitment (Phillis 2006, PMID: 16153659) [cell^2]")  # Table S4 (Kc_rec = 2.02e+07 cell^2)
    TGFbase <- fixed(0.08); label("Baseline TGFb level in breast tumor (Panis 2013, PMID: 23393376) [nmol/dm^3]")  # Table S4 (TGFbase = 8e-11 molarity)
    k_IFNg_sec <- fixed(4e-13); label("IFNg secretion rate by T helper cell (Borj 2017, doi:10.15171/ijbms.2017.05) [nmol/day*cell]")  # Table S4 (k_IFNg_sec = 4e-13 nanomole/cell/day)
    k_IFNg_deg <- fixed(11); label("IFNg degradation rate (Hofstra 1998, PMID: 9806748) [1/day]")  # Table S4 (k_IFNg_deg = 11 1/day)
    IFNg_50_ind <- fixed(0.00296); label("Half-Maximal IFNg level for PD-L1 induction (Shin 2017, PMID: 27903500) [nmol/dm^3]")  # Table S4 (IFNg_50_ind = 2.96e-12 molarity)
    r_carb <- fixed(1); label("Tumour to plasma concentration ratio of carboplatin (Konings 2009, PMID: 19096847) [1]")  # Table S4 (r_carb = 1 dimensionless)
    k_C_carb <- fixed(0.016728); label("Cancer cell killing rate by nab-paclitaxel (Wong 2018, PMID: 30323890) [1/day]")  # Table S4 (k_C_carb = 0.000697 1/h)
    IC50_carb <- fixed(318000); label("Half-maximal nab-paclitaxel concentration for cancer cell killing (Sadhukha 2014, PMID: 24831091) [nmol/dm^3]")  # Table S4 (IC50_carb = 0.000318 molarity)
    MW_carb <- fixed(0.37125); label("Molecular weight of nab-paclitaxel (calculated) [ug/nmol]")  # Table S4 (MW_carb = 371.25 gram/mole)
    Kc_carb <- fixed(80000000); label("Half-Maximal cancer cell number for cytotoxic drug diffusion (estimated based on Petersen 2022,PMID: 36159212) [cell]")  # Table S4 (Kc_carb = 8e+07 cell)
    V1_carb <- fixed(198); label("steady-state average concentration for Carboplatin when AUC=6 mg/mL/mi (calculated based on Jotte 2020,PMID: 32302702) [ug/dm^3]")  # Table S4 (V1_carb = 0.000198 gram/liter)

    # Three rows of Table S4 are deliberately not carried into ini() because they
    # are not referenced by any reaction or rule in Tables S5-S6:
    #   A_syn = 37.8 micrometer^2  -- duplicates the syn_T_C1 / syn_T_APC
    #     compartment capacities in Table S2, which are what the rate laws use.
    #   lagP = 3.6472 hour, durP = 0.1797 hour -- oral absorption lag time and
    #     zero-order duration for entinostat, a SimBiology dosing attribute
    #     rather than an ODE parameter. Entinostat is an inherited limb of the
    #     parent platform and is not dosed in any IMpower131 regimen, so both
    #     are inert here. Set them on the dosing record if entinostat is dosed.

    # Equation 2 (maximal cancer-cell capacity) was back-solved exactly from the
    # stored C_max = 862890782185.996 cell in Table S4:
    #   C_max = (4/3)*pi*(D_T_max/2)^3 * DEN_T_cell
    #         with D_T_max = 20 cm and DEN_T_cell = 2.06e8 cell/mL
    #         -> 8.62890782e11 cell, matching all 9 reported digits.
    # Neither constant enters the ODEs: Table S6 rule 51 makes C_max the
    # vasculature state V_T.K, so the tabulated C_max is a stale stored value.
  })

  model({
    # ---------- parameter aliases ----------
    # rxode2 parses a bare (theta + theta) sub-expression as a mu-reference block and
    # refuses the model. These mechanistic rate laws legitimately add two fixed
    # constants, so each such parameter is copied to a local first. Values and
    # provenance are unchanged.
    a_EC50_ArgI_Treg <- EC50_ArgI_Treg
    a_EC50_CCL2_rec <- EC50_CCL2_rec
    a_IC50_ArgI_CTL <- IC50_ArgI_CTL
    a_IC50_ENT_ArgI <- IC50_ENT_ArgI
    a_IC50_ENT_C <- IC50_ENT_C
    a_IC50_ENT_CCL2 <- IC50_ENT_CCL2
    a_IC50_ENT_NO <- IC50_ENT_NO
    a_IC50_NO_CTL <- IC50_NO_CTL
    a_IC50_carb <- IC50_carb
    a_IC50_nabp <- IC50_nabp
    a_IC50_nabp_vas <- IC50_nabp_vas
    a_IC50_vas <- IC50_vas
    a_IFNg_50_ind <- IFNg_50_ind
    a_IL2_50 <- IL2_50
    a_IL2_50_Treg <- IL2_50_Treg
    a_Kc_ENT <- Kc_ENT
    a_Kc_rec <- Kc_rec
    a_Kcl <- Kcl
    a_Kt <- Kt
    a_N0 <- N0
    a_TCR_p0_tot <- TCR_p0_tot
    a_TCR_p1_tot <- TCR_p1_tot
    a_TGF_50_ctl <- TGF_50_ctl
    a_TGF_50_reg <- TGF_50_reg
    a_TGFbase <- TGFbase
    a_c0 <- c0
    a_c50 <- c50
    a_cell <- cell
    a_k_C1_death <- k_C1_death
    a_k_M1p0_TCR_off <- k_M1p0_TCR_off
    a_k_M1p0_TCR_p <- k_M1p0_TCR_p
    a_k_M1p1_TCR_off <- k_M1p1_TCR_off
    a_k_M1p1_TCR_p <- k_M1p1_TCR_p
    a_k_T1_death <- k_T1_death
    a_p0_50 <- p0_50
    a_p1_50 <- p1_50
    a_phi_M1p0_TCR <- phi_M1p0_TCR
    a_phi_M1p1_TCR <- phi_M1p1_TCR

    # ---------- initial tumour geometry (Equation 2 / Table S6 rule 1) ----------
    vol_V_T0 <- (pi/6)*initial_tumour_diameter^3
    C1_0 <- vol_V_T0*Ve_T/vol_cell

    # ---------- algebraic definitions, in dependency order ----------
    m_V_C <- vol_V_C
    m_V_P <- vol_V_P
    m_V_LN <- vol_V_LN
    m_V_e <- vol_V_e
    m_A_e <- vol_A_e
    m_A_s <- vol_A_s
    m_syn_T_C1 <- vol_syn_T_C1
    m_syn_T_APC <- vol_syn_T_APC
    m_V_1 <- vol_V_1
    m_V_2 <- vol_V_2
    m_V_3 <- vol_V_3
    x_V_C_nT0 <- q_V_C_nT0
    x_V_C_T0 <- q_V_C_T0
    x_V_C_nT1 <- q_V_C_nT1
    x_V_C_T1 <- q_V_C_T1
    x_V_C_nivo <- q_V_C_nivo/m_V_C
    x_V_C_atezo <- q_V_C_atezo/m_V_C
    x_V_C_ipi <- q_V_C_ipi/m_V_C
    x_V_C_ENT <- q_V_C_ENT/m_V_C
    x_V_C_ENT_Buccal <- q_V_C_ENT_Buccal/m_V_C
    x_V_C_ENT_GI <- q_V_C_ENT_GI/m_V_C
    x_V_C_Dose2 <- q_V_C_Dose2/m_V_C
    x_V_C_Th <- q_V_C_Th
    x_V_P_nT0 <- q_V_P_nT0
    x_V_P_T0 <- q_V_P_T0
    x_V_P_nT1 <- q_V_P_nT1
    x_V_P_T1 <- q_V_P_T1
    x_V_P_nivo <- q_V_P_nivo/m_V_P
    x_V_P_atezo <- q_V_P_atezo/m_V_P
    x_V_P_ipi <- q_V_P_ipi/m_V_P
    x_V_P_Treg_CTLA4 <- q_V_P_Treg_CTLA4/NAVG_nmol
    x_V_P_Treg_CTLA4_ipi <- q_V_P_Treg_CTLA4_ipi/NAVG_nmol
    x_V_P_Treg_CTLA4_ipi_CTLA4 <- q_V_P_Treg_CTLA4_ipi_CTLA4/NAVG_nmol
    x_V_P_ENT <- q_V_P_ENT/m_V_P
    x_V_P_Th <- q_V_P_Th
    x_V_T_C_x <- q_V_T_C_x
    x_V_T_T1_exh <- q_V_T_T1_exh
    x_V_T_Th_exh <- q_V_T_Th_exh
    x_V_T_C1 <- q_V_T_C1
    x_V_T_C2 <- q_V_T_C2
    x_V_T_nT0 <- q_V_T_nT0
    x_V_T_T0 <- q_V_T_T0
    x_V_T_nT1 <- q_V_T_nT1
    x_V_T_T1 <- q_V_T_T1
    x_V_T_APC <- q_V_T_APC
    x_V_T_mAPC <- q_V_T_mAPC
    x_V_T_Treg_CTLA4 <- q_V_T_Treg_CTLA4/NAVG_nmol
    x_V_T_Treg_CTLA4_ipi <- q_V_T_Treg_CTLA4_ipi/NAVG_nmol
    x_V_T_Treg_CTLA4_ipi_CTLA4 <- q_V_T_Treg_CTLA4_ipi_CTLA4/NAVG_nmol
    x_V_T_MDSC <- q_V_T_MDSC
    x_V_T_ArgI <- q_V_T_ArgI
    x_V_T_K <- q_V_T_K
    x_V_T_Th <- q_V_T_Th
    x_V_LN_nT0 <- q_V_LN_nT0
    x_V_LN_aT0 <- q_V_LN_aT0
    x_V_LN_T0 <- q_V_LN_T0
    x_V_LN_IL2 <- q_V_LN_IL2/m_V_LN
    x_V_LN_nT1 <- q_V_LN_nT1
    x_V_LN_aT1 <- q_V_LN_aT1
    x_V_LN_T1 <- q_V_LN_T1
    x_V_LN_APC <- q_V_LN_APC
    x_V_LN_mAPC <- q_V_LN_mAPC
    x_V_LN_nivo <- q_V_LN_nivo/m_V_LN
    x_V_LN_atezo <- q_V_LN_atezo/m_V_LN
    x_V_LN_ipi <- q_V_LN_ipi/m_V_LN
    x_V_LN_ENT <- q_V_LN_ENT/m_V_LN
    x_V_LN_aTh <- q_V_LN_aTh
    x_V_LN_Th <- q_V_LN_Th
    x_V_e_P0 <- q_V_e_P0/m_V_e
    x_V_e_p0 <- q_V_e_p0/m_V_e
    x_V_e_P1 <- q_V_e_P1/m_V_e
    x_V_e_p1 <- q_V_e_p1/m_V_e
    x_A_e_M1 <- q_A_e_M1/NAVG_nmol/m_A_e
    x_A_e_M1p0 <- q_A_e_M1p0/NAVG_nmol/m_A_e
    x_A_e_M1p1 <- q_A_e_M1p1/NAVG_nmol/m_A_e
    x_A_s_M1 <- q_A_s_M1/NAVG_nmol/m_A_s
    x_A_s_M1p0 <- q_A_s_M1p0/NAVG_nmol/m_A_s
    x_A_s_M1p1 <- q_A_s_M1p1/NAVG_nmol/m_A_s
    x_syn_T_C1_PD1_PDL1 <- q_syn_T_C1_PD1_PDL1/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_PD1_PDL2 <- q_syn_T_C1_PD1_PDL2/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_PD1 <- q_syn_T_C1_PD1/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_PDL1 <- q_syn_T_C1_PDL1/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_PDL2 <- q_syn_T_C1_PDL2/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_PD1_nivo <- q_syn_T_C1_PD1_nivo/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_PD1_nivo_PD1 <- q_syn_T_C1_PD1_nivo_PD1/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_PDL1_atezo <- q_syn_T_C1_PDL1_atezo/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_PDL1_atezo_PDL1 <- q_syn_T_C1_PDL1_atezo_PDL1/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_TPDL1 <- q_syn_T_C1_TPDL1/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_TPDL1_atezo <- q_syn_T_C1_TPDL1_atezo/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_TPDL1_atezo_TPDL1 <- q_syn_T_C1_TPDL1_atezo_TPDL1/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_CD28_CD80 <- q_syn_T_C1_CD28_CD80/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_CD28_CD80_CD28 <- q_syn_T_C1_CD28_CD80_CD28/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_CD28_CD86 <- q_syn_T_C1_CD28_CD86/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_CD80_CTLA4 <- q_syn_T_C1_CD80_CTLA4/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_CD80_CTLA4_CD80 <- q_syn_T_C1_CD80_CTLA4_CD80/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_CTLA4_CD80_CTLA4 <- q_syn_T_C1_CTLA4_CD80_CTLA4/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_CD80_CTLA4_CD80_CTLA4 <- q_syn_T_C1_CD80_CTLA4_CD80_CTLA4/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_CD86_CTLA4 <- q_syn_T_C1_CD86_CTLA4/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_CD86_CTLA4_CD86 <- q_syn_T_C1_CD86_CTLA4_CD86/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_TPDL1_CD80 <- q_syn_T_C1_TPDL1_CD80/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_TPDL1_CD80_TPDL1 <- q_syn_T_C1_TPDL1_CD80_TPDL1/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_CD28 <- q_syn_T_C1_CD28/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_CTLA4 <- q_syn_T_C1_CTLA4/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_CD80 <- q_syn_T_C1_CD80/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_CD86 <- q_syn_T_C1_CD86/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_CTLA4_ipi <- q_syn_T_C1_CTLA4_ipi/NAVG_nmol/m_syn_T_C1
    x_syn_T_C1_CTLA4_ipi_CTLA4 <- q_syn_T_C1_CTLA4_ipi_CTLA4/NAVG_nmol/m_syn_T_C1
    x_syn_T_APC_PD1_PDL1 <- q_syn_T_APC_PD1_PDL1/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_PD1_PDL2 <- q_syn_T_APC_PD1_PDL2/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_PD1 <- q_syn_T_APC_PD1/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_PDL1 <- q_syn_T_APC_PDL1/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_PDL2 <- q_syn_T_APC_PDL2/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_PD1_nivo <- q_syn_T_APC_PD1_nivo/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_PD1_nivo_PD1 <- q_syn_T_APC_PD1_nivo_PD1/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_PDL1_atezo <- q_syn_T_APC_PDL1_atezo/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_PDL1_atezo_PDL1 <- q_syn_T_APC_PDL1_atezo_PDL1/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_TPDL1 <- q_syn_T_APC_TPDL1/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_TPDL1_atezo <- q_syn_T_APC_TPDL1_atezo/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_TPDL1_atezo_TPDL1 <- q_syn_T_APC_TPDL1_atezo_TPDL1/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_CD28_CD80 <- q_syn_T_APC_CD28_CD80/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_CD28_CD80_CD28 <- q_syn_T_APC_CD28_CD80_CD28/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_CD28_CD86 <- q_syn_T_APC_CD28_CD86/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_CD80_CTLA4 <- q_syn_T_APC_CD80_CTLA4/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_CD80_CTLA4_CD80 <- q_syn_T_APC_CD80_CTLA4_CD80/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_CTLA4_CD80_CTLA4 <- q_syn_T_APC_CTLA4_CD80_CTLA4/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_CD80_CTLA4_CD80_CTLA4 <- q_syn_T_APC_CD80_CTLA4_CD80_CTLA4/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_CD86_CTLA4 <- q_syn_T_APC_CD86_CTLA4/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_CD86_CTLA4_CD86 <- q_syn_T_APC_CD86_CTLA4_CD86/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_TPDL1_CD80 <- q_syn_T_APC_TPDL1_CD80/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_TPDL1_CD80_TPDL1 <- q_syn_T_APC_TPDL1_CD80_TPDL1/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_CD28 <- q_syn_T_APC_CD28/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_CTLA4 <- q_syn_T_APC_CTLA4/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_CD80 <- q_syn_T_APC_CD80/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_CD86 <- q_syn_T_APC_CD86/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_CTLA4_ipi <- q_syn_T_APC_CTLA4_ipi/NAVG_nmol/m_syn_T_APC
    x_syn_T_APC_CTLA4_ipi_CTLA4 <- q_syn_T_APC_CTLA4_ipi_CTLA4/NAVG_nmol/m_syn_T_APC
    x_V_1_NabP <- q_V_1_NabP/m_V_1
    x_V_2_NabP <- q_V_2_NabP/m_V_2
    x_V_3_NabP <- q_V_3_NabP/m_V_3
    C_total <- (((0*a_cell)+x_V_T_C1)+x_V_T_C2)   # Table S6 rule 2
    T_total <- ((((0*a_cell)+x_V_T_T0)+x_V_T_T1)+x_V_T_Th)   # Table S6 rule 3
    T_total_LN <- ((((0*a_cell)+x_V_LN_T0)+x_V_LN_T1)+x_V_LN_Th)   # Table S6 rule 4
    Tregs_ <- x_V_T_T0   # Table S6 rule 6
    H_APC <- ((n_sites_APC*x_V_LN_mAPC)/(((n_sites_APC*x_V_LN_mAPC)+(x_V_LN_nT0*n_T0_clones))+a_cell))   # Table S6 rule 9
    H_mAPC <- ((n_sites_APC*x_V_LN_mAPC)/(((n_sites_APC*x_V_LN_mAPC)+(x_V_LN_nT1*n_T1_clones))+a_cell))   # Table S6 rule 10
    H_APCh <- ((n_sites_APC*x_V_LN_mAPC)/(((n_sites_APC*x_V_LN_mAPC)+(x_V_LN_nT0*n_T1_clones))+a_cell))   # Table S6 rule 11
    pTCR_p0_MHC_tot <- (((a_k_M1p0_TCR_off/(a_k_M1p0_TCR_off+a_phi_M1p0_TCR))*((a_k_M1p0_TCR_p/(a_k_M1p0_TCR_off+a_k_M1p0_TCR_p))^N_M1p0_TCR))*(0.5*((((x_A_s_M1p0/n_T0_clones)+a_TCR_p0_tot)+(a_k_M1p0_TCR_off/k_M1p0_TCR_on))-(a_TCR_p0_tot*sqrt(((((((x_A_s_M1p0/n_T0_clones)+a_TCR_p0_tot)+(a_k_M1p0_TCR_off/k_M1p0_TCR_on))/a_TCR_p0_tot)^2)-(((4*x_A_s_M1p0)/n_T0_clones)/a_TCR_p0_tot)))))))   # Table S6 rule 12
    H_P0 <- (pTCR_p0_MHC_tot/(pTCR_p0_MHC_tot+a_p0_50))   # Table S6 rule 13
    pTCR_p1_MHC_tot <- (((a_k_M1p1_TCR_off/(a_k_M1p1_TCR_off+a_phi_M1p1_TCR))*((a_k_M1p1_TCR_p/(a_k_M1p1_TCR_off+a_k_M1p1_TCR_p))^N_M1p1_TCR))*(0.5*((((x_A_s_M1p1/n_T1_clones)+a_TCR_p1_tot)+(a_k_M1p1_TCR_off/k_M1p1_TCR_on))-(a_TCR_p1_tot*sqrt(((((((x_A_s_M1p1/n_T1_clones)+a_TCR_p1_tot)+(a_k_M1p1_TCR_off/k_M1p1_TCR_on))/a_TCR_p1_tot)^2)-(((4*x_A_s_M1p1)/n_T1_clones)/a_TCR_p1_tot)))))))   # Table S6 rule 14
    H_P1 <- (pTCR_p1_MHC_tot/(pTCR_p1_MHC_tot+a_p1_50))   # Table S6 rule 15
    x_syn_T_C1_PDL1_total <- (((x_syn_T_C1_PDL1+x_syn_T_C1_PD1_PDL1)+x_syn_T_C1_PDL1_atezo)+(2*x_syn_T_C1_PDL1_atezo_PDL1))   # Table S6 rule 22
    x_syn_T_C1_PDL2_total <- (x_syn_T_C1_PD1_PDL2+x_syn_T_C1_PDL2)   # Table S6 rule 23
    H_PD1_C1 <- ((((x_syn_T_C1_PD1_PDL1+x_syn_T_C1_PD1_PDL2)/PD1_50)^n_PD1)/((((x_syn_T_C1_PD1_PDL1+x_syn_T_C1_PD1_PDL2)/PD1_50)^n_PD1)+1))   # Table S6 rule 26
    H_CD28_C1 <- (((((x_syn_T_C1_CD28_CD80+x_syn_T_C1_CD28_CD86)+(2*x_syn_T_C1_CD28_CD80_CD28))/CD28_CD8X_50)^n_CD28_CD8X)/(((((x_syn_T_C1_CD28_CD80+x_syn_T_C1_CD28_CD86)+(2*x_syn_T_C1_CD28_CD80_CD28))/CD28_CD8X_50)^n_CD28_CD8X)+1))   # Table S6 rule 27
    x_syn_T_APC_PDL1_total <- (((x_syn_T_APC_PDL1+x_syn_T_APC_PD1_PDL1)+x_syn_T_APC_PDL1_atezo)+(2*x_syn_T_APC_PDL1_atezo_PDL1))   # Table S6 rule 34
    x_syn_T_APC_PDL2_total <- (x_syn_T_APC_PD1_PDL2+x_syn_T_APC_PDL2)   # Table S6 rule 35
    H_PD1_APC <- ((((x_syn_T_APC_PD1_PDL1+x_syn_T_APC_PD1_PDL2)/PD1_50)^n_PD1)/((((x_syn_T_APC_PD1_PDL1+x_syn_T_APC_PD1_PDL2)/PD1_50)^n_PD1)+1))   # Table S6 rule 38
    H_CD28_APC <- (((((x_syn_T_APC_CD28_CD80+x_syn_T_APC_CD28_CD86)+(2*x_syn_T_APC_CD28_CD80_CD28))/CD28_CD8X_50)^n_CD28_CD8X)/(((((x_syn_T_APC_CD28_CD80+x_syn_T_APC_CD28_CD86)+(2*x_syn_T_APC_CD28_CD80_CD28))/CD28_CD8X_50)^n_CD28_CD8X)+1))   # Table S6 rule 39
    H_Treg_T <- ((((x_V_T_Treg_CTLA4_ipi+(2*x_V_T_Treg_CTLA4_ipi_CTLA4))/Treg_CTLA4_50)^n_Treg_CTLA4)/((((x_V_T_Treg_CTLA4_ipi+(2*x_V_T_Treg_CTLA4_ipi_CTLA4))/Treg_CTLA4_50)^n_Treg_CTLA4)+1))   # Table S6 rule 42
    H_Treg_P <- ((((x_V_P_Treg_CTLA4_ipi+(2*x_V_P_Treg_CTLA4_ipi_CTLA4))/Treg_CTLA4_50)^n_Treg_CTLA4)/((((x_V_P_Treg_CTLA4_ipi+(2*x_V_P_Treg_CTLA4_ipi_CTLA4))/Treg_CTLA4_50)^n_Treg_CTLA4)+1))   # Table S6 rule 43
    x_V_T_NabP <- ((r_nabp*x_V_1_NabP)/MW_nabp)   # Table S6 rule 46
    k_C1_therapy <- (((k_C_nabp*(x_V_T_NabP/(x_V_T_NabP+a_IC50_nabp)))*min(C_total, Kc_cyt))/C_total)   # Table S6 rule 47
    k_C2_therapy <- (((k_C_nabp*(x_V_T_NabP/(x_V_T_NabP+(a_IC50_nabp*r_resist))))*min(C_total, Kc_cyt))/C_total)   # Table S6 rule 48
    k_C2_growth <- k_C1_growth   # Table S6 rule 49
    k_C2_death <- a_k_C1_death   # Table S6 rule 50
    C_max <- x_V_T_K   # Table S6 rule 51
    N_aTh <- ((a_N0+(N_costim*H_CD28_APC))+((N_IL2_CD4*x_V_LN_IL2)/(a_IL2_50+x_V_LN_IL2)))   # Table S6 rule 54
    x_V_T_carb <- ((r_carb*V1_carb)/MW_carb)   # Table S6 rule 55
    vol_V_T <- ((((x_V_T_C_x+C_total)*vol_cell)+((((x_V_T_T1_exh+x_V_T_Th_exh)+T_total)+x_V_T_Th)*vol_Tcell))/Ve_T)   # Table S6 rule 1
    N_aT <- ((a_N0+(N_costim*H_CD28_APC))+((N_IL2_CD8*x_V_LN_IL2)/(a_IL2_50+x_V_LN_IL2)))   # Table S6 rule 7
    N_aT0 <- ((a_N0+(N_costim*H_CD28_APC))+((N_IL2_CD4*x_V_LN_IL2)/(a_IL2_50_Treg+x_V_LN_IL2)))   # Table S6 rule 8
    x_V_T_c <- q_V_T_c/vol_V_T
    x_V_T_P0 <- q_V_T_P0/vol_V_T
    x_V_T_P1 <- q_V_T_P1/vol_V_T
    x_V_T_nivo <- q_V_T_nivo/vol_V_T
    x_V_T_atezo <- q_V_T_atezo/vol_V_T
    x_V_T_ipi <- q_V_T_ipi/vol_V_T
    x_V_T_CCL2 <- q_V_T_CCL2/vol_V_T
    x_V_T_NO <- q_V_T_NO/vol_V_T
    x_V_T_ENT <- q_V_T_ENT/vol_V_T
    x_V_T_c_vas <- q_V_T_c_vas/vol_V_T
    x_V_T_TGF <- q_V_T_TGF/vol_V_T
    x_V_T_IFNg <- q_V_T_IFNg/vol_V_T
    H_MDSC_C1 <- (1-((a_IC50_NO_CTL/(a_IC50_NO_CTL+x_V_T_NO))*(1-((x_V_T_ArgI/(a_IC50_ArgI_CTL+x_V_T_ArgI))*(a_IC50_ENT_ArgI/(x_V_T_ENT+a_IC50_ENT_ArgI))))))   # Table S6 rule 44
    H_ENT_C1 <- (x_V_T_ENT/(x_V_T_ENT+a_IC50_ENT_C))   # Table S6 rule 45
    H_TGF_reg <- (x_V_T_TGF/(x_V_T_TGF+a_TGF_50_reg))   # Table S6 rule 52
    H_TGF_CTL <- (x_V_T_TGF/(x_V_T_TGF+a_TGF_50_ctl))   # Table S6 rule 53
    R_Tcell <- ((((((0*a_cell)/day)+((a_k_T1_death+k_C1_therapy)*x_V_T_C1))+((((((((k_C_T1*x_V_T_C1)*x_V_T_T1)/(((1.2*x_V_T_C1)+x_V_T_T1)+a_cell))*x_V_T_T1)/((x_V_T_T1+(11*Tregs_))+a_cell))*(1-H_TGF_CTL))*(1-H_PD1_C1))*(1-H_MDSC_C1)))+((a_k_T1_death+k_C2_therapy)*x_V_T_C2))+((((((((k_C_T1*x_V_T_C2)*x_V_T_T1)/(((1.2*x_V_T_C2)+x_V_T_T1)+a_cell))*x_V_T_T1)/((x_V_T_T1+(11*Tregs_))+a_cell))*(1-H_TGF_CTL))*(1-H_PD1_C1))*(1-H_MDSC_C1)))   # Table S6 rule 5

    # ---------- reaction fluxes (Supplementary Table S5) ----------
    v1 <- (k_cell_clear*x_V_T_C_x)   # R1 V_T.C_x -> null: Clearance of dead cancer cells from tumor
    v2 <- (k_cell_clear*x_V_T_T1_exh)   # R2 V_T.T1_exh -> null: Clearance of exhausted CD8 T cells from tumor
    v3 <- (k_cell_clear*x_V_T_Th_exh)   # R3 V_T.Th_exh -> null: Clearance of exhausted CD4 T cells from tumor
    v4 <- (((k_C1_growth*x_V_T_C1)*log(max((C_max/(C_total+a_cell)), 1)))*(1-H_ENT_C1))   # R4 null -> V_T.C1: Cancer cell growth
    v5 <- (a_k_C1_death*x_V_T_C1)   # R5 V_T.C1 -> V_T.C_x: Cancer cell death
    v6 <- (((k_C2_growth*x_V_T_C2)*log(max((C_max/(C_total+a_cell)), 1)))*(1-H_ENT_C1))   # R6 null -> V_T.C2: Cancer cell growth
    v7 <- (k_C2_death*x_V_T_C2)   # R7 V_T.C2 -> V_T.C_x: Cancer cell death
    v8 <- (Q_nT0_thym/div_T0)   # R8 null -> V_C.nT0: Thymic output of naive T cell to blood
    v9 <- (((k_nT0_pro/div_T0)*x_V_P_nT0)/((K_nT0_pro/div_T0)+x_V_P_nT0))   # R9 null -> V_P.nT0: Naive T cell proliferation in the peripheral compartment
    v10 <- (((k_nT0_pro/div_T0)*x_V_LN_nT0)/((K_nT0_pro/div_T0)+x_V_LN_nT0))   # R10 null -> V_LN.nT0: Naive T cell proliferation in the TDLN compartment
    v11 <- (k_nT0_death*x_V_P_nT0)   # R11 V_P.nT0 -> null: Naive T cell death in the peripheral compartment
    v12 <- (k_nT0_death*x_V_C_nT0)   # R12 V_C.nT0 -> null: Naive T cell death in the central compartment
    v13 <- (k_nT0_death*x_V_LN_nT0)   # R13 V_LN.nT0 -> null: Naive T cell death in the TDLN compartment
    v14 <- (k_nT0_death*x_V_T_nT0)   # R14 V_T.nT0 -> null: Naive T cell death in the tumor compartment
    v15 <- (q_nT0_P_in*x_V_C_nT0)   # R15 V_C.nT0 -> V_P.nT0: Naive T cell entry into the peripheral compartment
    v16 <- (q_nT0_P_out*x_V_P_nT0)   # R16 V_P.nT0 -> V_C.nT0: Naive T cell exit from the peripheral compartment
    v17 <- (q_nT0_LN_in*x_V_C_nT0)   # R17 V_C.nT0 -> V_LN.nT0: Naive T cell entry into the lymph node
    v18 <- (q_nT0_LN_out*x_V_LN_nT0)   # R18 V_LN.nT0 -> V_C.nT0: Naive T cell exit from the lymph node
    v19 <- (((k_T0_act*H_APC)*H_P0)*x_V_LN_nT0)   # R19 V_LN.nT0 -> null: Naive T cell activation
    v20 <- ((((k_T0_act*H_APC)*H_P0)*x_V_LN_nT0)*n_T0_clones)   # R20 null -> V_LN.aT0: Naive T cell activation
    v21 <- ((k_T0_pro/N_aT0)*x_V_LN_aT0)   # R21 V_LN.aT0 -> null: aT0 cell proliferation
    v22 <- (((k_T0_pro/N_aT0)*(2^N_aT0))*x_V_LN_aT0)   # R22 null -> V_LN.T0: aT0 cell proliferation
    v23 <- (k_T0_death*x_V_C_T0)   # R23 V_C.T0 -> null: T cell death in the central compartment
    v24 <- (k_T0_death*x_V_P_T0)   # R24 V_P.T0 -> null: T cell death in the peripheral compartment
    v25 <- (k_T0_death*x_V_LN_T0)   # R25 V_LN.T0 -> null: T cell death in the lymph node compartment
    v26 <- (k_T0_death*x_V_T_T0)   # R26 V_T.T0 -> null: T cell death in the tumor compartment
    v27 <- (q_T0_P_in*x_V_C_T0)   # R27 V_C.T0 -> V_P.T0: T cell transport into the peripheral compartment
    v28 <- (q_T0_P_out*x_V_P_T0)   # R28 V_P.T0 -> V_C.T0: T cell transport out of the peripheral compartment
    v29 <- (((q_T0_T_in*vol_V_T)*x_V_C_T0)*((C_total^2)/((C_total^2)+a_Kc_rec)))   # R29 V_C.T0 -> V_T.T0: T cell transport into the tumor compartment
    v30 <- (q_T0_LN_out*x_V_LN_T0)   # R30 V_LN.T0 -> V_C.T0: T cell transport out of the lymph node compartment
    v31 <- ((k_IL2_deg*x_V_LN_IL2)*m_V_LN)   # R31 V_LN.IL2 -> null: IL2 degradation
    v32 <- (((k_IL2_cons*x_V_LN_T1)*x_V_LN_IL2)/(a_IL2_50+x_V_LN_IL2))   # R32 V_LN.IL2 -> null: IL2 consumption by T cells
    v33 <- (((k_IL2_cons*x_V_LN_T0)*x_V_LN_IL2)/(a_IL2_50_Treg+x_V_LN_IL2))   # R33 V_LN.IL2 -> null: IL2 consumption by Tregs
    v34 <- (Q_nT1_thym/div_T1)   # R34 null -> V_C.nT1: Thymic output of naive T cell to blood
    v35 <- (((k_nT1_pro/div_T1)*x_V_P_nT1)/((K_nT1_pro/div_T1)+x_V_P_nT1))   # R35 null -> V_P.nT1: Naive T cell proliferation in the peripheral compartment
    v36 <- (((k_nT1_pro/div_T1)*x_V_LN_nT1)/((K_nT1_pro/div_T1)+x_V_LN_nT1))   # R36 null -> V_LN.nT1: Naive T cell proliferation in the TDLN compartment
    v37 <- (k_nT1_death*x_V_P_nT1)   # R37 V_P.nT1 -> null: Naive T cell death in the peripheral compartment
    v38 <- (k_nT1_death*x_V_C_nT1)   # R38 V_C.nT1 -> null: Naive T cell death in the central compartment
    v39 <- (k_nT1_death*x_V_LN_nT1)   # R39 V_LN.nT1 -> null: Naive T cell death in the TDLN compartment
    v40 <- (k_nT1_death*x_V_T_nT1)   # R40 V_T.nT1 -> null: Naive T cell death in the tumor compartment
    v41 <- (q_nT1_P_in*x_V_C_nT1)   # R41 V_C.nT1 -> V_P.nT1: Naive T cell entry into the peripheral compartment
    v42 <- (q_nT1_P_out*x_V_P_nT1)   # R42 V_P.nT1 -> V_C.nT1: Naive T cell exit from the peripheral compartment
    v43 <- (q_nT1_LN_in*x_V_C_nT1)   # R43 V_C.nT1 -> V_LN.nT1: Naive T cell entry into the lymph node
    v44 <- (q_nT1_LN_out*x_V_LN_nT1)   # R44 V_LN.nT1 -> V_C.nT1: Naive T cell exit from the lymph node
    v45 <- (((k_T1_act*H_mAPC)*H_P1)*x_V_LN_nT1)   # R45 V_LN.nT1 -> null: Naive T cell activation
    v46 <- ((((k_T1_act*H_mAPC)*H_P1)*x_V_LN_nT1)*n_T1_clones)   # R46 null -> V_LN.aT1: Naive T cell activation
    v47 <- ((k_T1_pro/N_aT)*x_V_LN_aT1)   # R47 V_LN.aT1 -> null: aT1 cell proliferation
    v48 <- (((k_T1_pro/N_aT)*(2^N_aT))*x_V_LN_aT1)   # R48 null -> V_LN.T1: aT1 cell proliferation
    v49 <- (a_k_T1_death*x_V_C_T1)   # R49 V_C.T1 -> null: T cell death in the central compartment
    v50 <- (a_k_T1_death*x_V_P_T1)   # R50 V_P.T1 -> null: T cell death in the peripheral compartment
    v51 <- (a_k_T1_death*x_V_LN_T1)   # R51 V_LN.T1 -> null: T cell death in the lymph node compartment
    v52 <- (a_k_T1_death*x_V_T_T1)   # R52 V_T.T1 -> V_T.T1_exh: T cell death in the tumor compartment
    v53 <- (((k_Treg*x_V_T_T1)*Tregs_)/((x_V_T_T1+Tregs_)+a_cell))   # R53 V_T.T1 -> V_T.T1_exh: T cell death from Tregs
    v54 <- ((((k_T1*x_V_T_T1)*C_total)/((C_total+x_V_T_T1)+a_cell))*H_PD1_C1)   # R54 V_T.T1 -> V_T.T1_exh: T cell death from cancer
    v55 <- (q_T1_P_in*x_V_C_T1)   # R55 V_C.T1 -> V_P.T1: T cell transport into the peripheral compartment
    v56 <- (q_T1_P_out*x_V_P_T1)   # R56 V_P.T1 -> V_C.T1: T cell transport out of the peripheral compartment
    v57 <- (((q_T1_T_in*vol_V_T)*x_V_C_T1)*((C_total^2)/((C_total^2)+a_Kc_rec)))   # R57 V_C.T1 -> V_T.T1: T cell transport into the tumor compartment
    v58 <- (q_T1_LN_out*x_V_LN_T1)   # R58 V_LN.T1 -> V_C.T1: T cell transport out of the lymph node compartment
    v59 <- (k_IL2_sec*x_V_LN_aT1)   # R59 null -> V_LN.IL2: IL2 secretion from activated T cells
    v60 <- ((((((((k_C_T1*x_V_T_C1)*x_V_T_T1)/(((1.2*x_V_T_C1)+x_V_T_T1)+a_cell))*x_V_T_T1)/((x_V_T_T1+(11*Tregs_))+a_cell))*(1-H_TGF_CTL))*(1-H_PD1_C1))*(1-H_MDSC_C1))   # R60 V_T.C1 -> V_T.C_x: Cancer cell killing by T cells
    v61 <- ((((((((k_C_T1*x_V_T_C2)*x_V_T_T1)/(((1.2*x_V_T_C2)+x_V_T_T1)+a_cell))*x_V_T_T1)/((x_V_T_T1+(11*Tregs_))+a_cell))*(1-H_TGF_CTL))*(1-H_PD1_C1))*(1-H_MDSC_C1))   # R61 V_T.C2 -> V_T.C_x: Cancer cell killing by T cells
    v62 <- (k_APC_death*((APC0_T*vol_V_T)-x_V_T_APC))   # R62 null -> V_T.APC: APC recruitment/death in the tumour
    v63 <- (k_APC_death*((APC0_LN*m_V_LN)-x_V_LN_APC))   # R63 null -> V_LN.APC: APC recruitment/death in LN
    v64 <- (((k_APC_mat*x_V_T_c)/(x_V_T_c+a_c50))*x_V_T_APC)   # R64 V_T.APC -> V_T.mAPC: APC maturation in the tumour
    v65 <- (k_APC_mig*x_V_T_mAPC)   # R65 V_T.mAPC -> V_LN.mAPC: APC migration to the lymph node
    v66 <- (k_mAPC_death*x_V_T_mAPC)   # R66 V_T.mAPC -> null: mAPC death in the tumour
    v67 <- (k_mAPC_death*x_V_LN_mAPC)   # R67 V_LN.mAPC -> null: mAPC death in the lymph node
    v68 <- (k_c*(a_c0-x_V_T_c))   # R68 null -> V_T.c: Baseline cytokine secretion/degradation
    v69 <- (R_Tcell*DAMPs)   # R69 null -> V_T.c: Cytokine release in response to the tumour
    v70 <- (((kout*x_A_e_M1)*m_A_e)-((kin*x_A_s_M1)*m_A_s))   # R70 A_e.M1 -> A_s.M1: MHC translocation
    v71 <- (n_T0_clones*(((((P0_C1*(a_k_C1_death+k_C1_therapy))*x_V_T_C1)+(((((((((P0_C1*k_C_T1)*x_V_T_C1)*x_V_T_T1)/(((1.2*x_V_T_C1)+x_V_T_T1)+a_cell))*x_V_T_T1)/((x_V_T_T1+(11*Tregs_))+a_cell))*(1-H_TGF_CTL))*(1-H_PD1_C1))*(1-H_MDSC_C1)))+((P0_C2*(k_C2_death+k_C2_therapy))*x_V_T_C2))+(((((((((P0_C2*k_C_T1)*x_V_T_C2)*x_V_T_T1)/(((1.2*x_V_T_C2)+x_V_T_T1)+a_cell))*x_V_T_T1)/((x_V_T_T1+(11*Tregs_))+a_cell))*(1-H_TGF_CTL))*(1-H_PD1_C1))*(1-H_MDSC_C1))))   # R71 null -> V_T.P0: Antigen deposition from dying cancer cells
    v72 <- ((k_xP0_deg*x_V_T_P0)*vol_V_T)   # R72 V_T.P0 -> null: Free antigen degradation
    v73 <- (((k_P0_up*x_V_T_APC)*x_V_T_P0)*vol_V_T)   # R73 V_T.P0 -> null: Antigen uptake by mature antigen presenting cells
    v74 <- (((k_P0_up*a_cell)*x_V_T_P0)*m_V_e)   # R74 null -> V_e.P0: Antigen uptake by mature antigen presenting cells
    v75 <- ((k_P0_deg*x_V_e_P0)*m_V_e)   # R75 V_e.P0 -> V_e.p0: Antigen degradation in APC endosomes
    v76 <- ((k_p0_deg*x_V_e_p0)*m_V_e)   # R76 V_e.p0 -> null: Epitope degradation in APC endosomes
    v77 <- (((k_P0_on*x_V_e_p0)*x_A_e_M1)*m_A_e)   # R77 V_e.p0 + A_e.M1 -> A_e.M1p0: Antigen-MHC binding in endosome
    v78 <- (((k_P0_d1*k_P0_on)*x_A_e_M1p0)*m_A_e)   # R78 A_e.M1p0 -> V_e.p0 + A_e.M1: Antigen-MHC unbinding in endosome
    v79 <- (((k_P0_d1*k_P0_on)*x_A_s_M1p0)*m_A_s)   # R79 A_s.M1p0 -> A_s.M1: Antigen-MHC unbinding on APC surface
    v80 <- ((kout*x_A_e_M1p0)*m_A_e)   # R80 A_e.M1p0 -> A_s.M1p0: Antigen-MHC translocation
    v81 <- (n_T1_clones*(((((P1_C1*(a_k_C1_death+k_C1_therapy))*x_V_T_C1)+(((((((((P1_C1*k_C_T1)*x_V_T_C1)*x_V_T_T1)/(((1.2*x_V_T_C1)+x_V_T_T1)+a_cell))*x_V_T_T1)/((x_V_T_T1+(11*Tregs_))+a_cell))*(1-H_TGF_CTL))*(1-H_PD1_C1))*(1-H_MDSC_C1)))+((P1_C2*(k_C2_death+k_C2_therapy))*x_V_T_C2))+(((((((((P1_C2*k_C_T1)*x_V_T_C2)*x_V_T_T1)/(((1.2*x_V_T_C2)+x_V_T_T1)+a_cell))*x_V_T_T1)/((x_V_T_T1+(11*Tregs_))+a_cell))*(1-H_TGF_CTL))*(1-H_PD1_C1))*(1-H_MDSC_C1))))   # R81 null -> V_T.P1: Antigen deposition from dying cancer cells
    v82 <- ((k_xP1_deg*x_V_T_P1)*vol_V_T)   # R82 V_T.P1 -> null: Free antigen degradation
    v83 <- (((k_P1_up*x_V_T_APC)*x_V_T_P1)*vol_V_T)   # R83 V_T.P1 -> null: Antigen uptake by mature antigen presenting cells
    v84 <- (((k_P1_up*a_cell)*x_V_T_P1)*m_V_e)   # R84 null -> V_e.P1: Antigen uptake by mature antigen presenting cells
    v85 <- ((k_P1_deg*x_V_e_P1)*m_V_e)   # R85 V_e.P1 -> V_e.p1: Antigen degradation in APC endosomes
    v86 <- ((k_p1_deg*x_V_e_p1)*m_V_e)   # R86 V_e.p1 -> null: Epitope degradation in APC endosomes
    v87 <- (((k_P1_on*x_V_e_p1)*x_A_e_M1)*m_A_e)   # R87 V_e.p1 + A_e.M1 -> A_e.M1p1: Antigen-MHC binding in endosome
    v88 <- (((k_P1_d1*k_P1_on)*x_A_e_M1p1)*m_A_e)   # R88 A_e.M1p1 -> V_e.p1 + A_e.M1: Antigen-MHC unbinding in endosome
    v89 <- (((k_P1_d1*k_P1_on)*x_A_s_M1p1)*m_A_s)   # R89 A_s.M1p1 -> A_s.M1: Antigen-MHC unbinding on APC surface
    v90 <- ((kout*x_A_e_M1p1)*m_A_e)   # R90 A_e.M1p1 -> A_s.M1p1: Antigen-MHC translocation
    v91 <- (q_P_nivo*((x_V_C_nivo/gamma_C_nivo)-(x_V_P_nivo/gamma_P_nivo)))   # R91 V_C.nivo <-> V_P.nivo: nivo diffusive transport to peripheral compartment
    v92 <- (q_T_nivo*((x_V_C_nivo/gamma_C_nivo)-(x_V_T_nivo/gamma_T_nivo)))   # R92 V_C.nivo <-> V_T.nivo: nivo diffusive transport to tumor compartment
    v93 <- (q_LN_nivo*((x_V_C_nivo/gamma_C_nivo)-(x_V_LN_nivo/gamma_LN_nivo)))   # R93 V_C.nivo <-> V_LN.nivo: nivo diffusive transport to lymph node compartment
    v94 <- ((q_LD_nivo*x_V_T_nivo)/gamma_T_nivo)   # R94 V_T.nivo -> V_LN.nivo: nivo convective transport from tumor to lymph node
    v95 <- ((q_LD_nivo*x_V_LN_nivo)/gamma_LN_nivo)   # R95 V_LN.nivo -> V_C.nivo: nivo convective transport from lymph node to central
    v96 <- (k_cl_nivo*x_V_C_nivo)   # R96 V_C.nivo -> null: Drug clearance from central compartment
    v97 <- (q_P_atezo*((x_V_C_atezo/gamma_C_atezo)-(x_V_P_atezo/gamma_P_atezo)))   # R97 V_C.atezo <-> V_P.atezo: atezo diffusive transport to peripheral compartment
    v98 <- (q_T_atezo*((x_V_C_atezo/gamma_C_atezo)-(x_V_T_atezo/gamma_T_atezo)))   # R98 V_C.atezo <-> V_T.atezo: atezo diffusive transport to tumor compartment
    v99 <- (q_LN_atezo*((x_V_C_atezo/gamma_C_atezo)-(x_V_LN_atezo/gamma_LN_atezo)))   # R99 V_C.atezo <-> V_LN.atezo: atezo diffusive transport to lymph node compartment
    v100 <- ((q_LD_atezo*x_V_T_atezo)/gamma_T_atezo)   # R100 V_T.atezo -> V_LN.atezo: atezo convective transport from tumor to lymph node
    v101 <- ((q_LD_atezo*x_V_LN_atezo)/gamma_LN_atezo)   # R101 V_LN.atezo -> V_C.atezo: atezo convective transport from lymph node to central
    v102 <- (k_cl_atezo*x_V_C_atezo)   # R102 V_C.atezo -> null: Drug clearance from central compartment
    v103 <- (q_P_ipi*((x_V_C_ipi/gamma_C_ipi)-(x_V_P_ipi/gamma_P_ipi)))   # R103 V_C.ipi <-> V_P.ipi: ipi diffusive transport to peripheral compartment
    v104 <- (q_T_ipi*((x_V_C_ipi/gamma_C_ipi)-(x_V_T_ipi/gamma_T_ipi)))   # R104 V_C.ipi <-> V_T.ipi: ipi diffusive transport to tumor compartment
    v105 <- (q_LN_ipi*((x_V_C_ipi/gamma_C_ipi)-(x_V_LN_ipi/gamma_LN_ipi)))   # R105 V_C.ipi <-> V_LN.ipi: ipi diffusive transport to lymph node compartment
    v106 <- ((q_LD_ipi*x_V_T_ipi)/gamma_T_ipi)   # R106 V_T.ipi -> V_LN.ipi: ipi convective transport from tumor to lymph node
    v107 <- ((q_LD_ipi*x_V_LN_ipi)/gamma_LN_ipi)   # R107 V_LN.ipi -> V_C.ipi: ipi convective transport from lymph node to central
    v108 <- (k_cl_ipi*x_V_C_ipi)   # R108 V_C.ipi -> null: Drug clearance from central compartment
    v109 <- (((k_out_PDL1*x_V_T_IFNg)/(x_V_T_IFNg+a_IFNg_50_ind))*(1-(x_syn_T_C1_PDL1_total/((C1_PDL1_base*r_PDL1_IFNg)/A_cell))))   # R109 null -> syn_T_C1.PDL1: Translocation of PDL1 between cell surface and cytoplasm
    v110 <- ((((k_out_PDL1*r_PDL2C1)*x_V_T_IFNg)/(x_V_T_IFNg+a_IFNg_50_ind))*(1-(x_syn_T_C1_PDL2_total/(((C1_PDL1_base*r_PDL1_IFNg)/A_cell)*r_PDL2C1))))   # R110 null -> syn_T_C1.PDL2: Translocation of PDL2 between cell surface and cytoplasm
    v111 <- (k_in_PDL1*((C1_PDL1_base/A_cell)-x_syn_T_C1_PDL1_total))   # R111 null -> syn_T_C1.PDL1: Translocation of PDL1 between cell surface and cytoplasm
    v112 <- (k_in_PDL1*(((C1_PDL1_base/A_cell)*r_PDL2C1)-x_syn_T_C1_PDL2_total))   # R112 null -> syn_T_C1.PDL2: Translocation of PDL2 between cell surface and cytoplasm
    v113 <- (((kon_PD1_PDL1*x_syn_T_C1_PD1)*x_syn_T_C1_PDL1)-(koff_PD1_PDL1*x_syn_T_C1_PD1_PDL1))   # R113 syn_T_C1.PD1 + syn_T_C1.PDL1 <-> syn_T_C1.PD1_PDL1: binding and unbinding of PD1 PDL1 in synapse
    v114 <- (((kon_PD1_PDL2*x_syn_T_C1_PD1)*x_syn_T_C1_PDL2)-(koff_PD1_PDL2*x_syn_T_C1_PD1_PDL2))   # R114 syn_T_C1.PD1 + syn_T_C1.PDL2 <-> syn_T_C1.PD1_PDL2: binding and unbinding of PD1 PDL2 in synapse
    v115 <- (((2*kon_PD1_nivo)*((x_syn_T_C1_PD1*x_V_T_nivo)/gamma_T_nivo))-(koff_PD1_nivo*x_syn_T_C1_PD1_nivo))   # R115 syn_T_C1.PD1 <-> syn_T_C1.PD1_nivo: binding and unbinding of PD1 to Nivo on T surface in synapse
    v116 <- (((Chi_PD1_nivo*kon_PD1_nivo)*(x_syn_T_C1_PD1*x_syn_T_C1_PD1_nivo))-((2*koff_PD1_nivo)*x_syn_T_C1_PD1_nivo_PD1))   # R116 syn_T_C1.PD1_nivo + syn_T_C1.PD1 <-> syn_T_C1.PD1_nivo_PD1: binding and unbinding of PD1 to PD1-Nivo on T surface in synapse
    v117 <- (((2*kon_PDL1_atezo)*((x_syn_T_C1_PDL1*x_V_T_atezo)/gamma_T_atezo))-(koff_PDL1_atezo*x_syn_T_C1_PDL1_atezo))   # R117 syn_T_C1.PDL1 <-> syn_T_C1.PDL1_atezo: binding and unbinding of PDL1 to atezo on C1 surface in synapse
    v118 <- (((Chi_PDL1_atezo*kon_PDL1_atezo)*(x_syn_T_C1_PDL1*x_syn_T_C1_PDL1_atezo))-((2*koff_PDL1_atezo)*x_syn_T_C1_PDL1_atezo_PDL1))   # R118 syn_T_C1.PDL1_atezo + syn_T_C1.PDL1 <-> syn_T_C1.PDL1_atezo_PDL1: binding and unbinding of PDL1 to PDL1-atezo on C1 surface in synapse
    v119 <- ((((2*kon_CD28_CD80)*x_syn_T_C1_CD28)*x_syn_T_C1_CD80)-(koff_CD28_CD80*x_syn_T_C1_CD28_CD80))   # R119 syn_T_C1.CD28 + syn_T_C1.CD80 <-> syn_T_C1.CD28_CD80: binding and unbinding of CD28 and CD80 in synapse
    v120 <- (((kon_CD28_CD80*x_syn_T_C1_CD28)*x_syn_T_C1_CD28_CD80)-((2*koff_CD28_CD80)*x_syn_T_C1_CD28_CD80_CD28))   # R120 syn_T_C1.CD28_CD80 + syn_T_C1.CD28 <-> syn_T_C1.CD28_CD80_CD28: binding and unbinding of CD28-CD80 and CD28 in synapse
    v121 <- (((kon_CD28_CD86*x_syn_T_C1_CD28)*x_syn_T_C1_CD86)-(koff_CD28_CD86*x_syn_T_C1_CD28_CD86))   # R121 syn_T_C1.CD28 + syn_T_C1.CD86 <-> syn_T_C1.CD28_CD86: binding and unbinding of CD28 and CD86 in synapse
    v122 <- ((((4*kon_CTLA4_CD80)*x_syn_T_C1_CTLA4)*x_syn_T_C1_CD80)-(koff_CTLA4_CD80*x_syn_T_C1_CD80_CTLA4))   # R122 syn_T_C1.CTLA4 + syn_T_C1.CD80 <-> syn_T_C1.CD80_CTLA4: binding and unbinding of CTLA4 and CD80 in synapse
    v123 <- (((kon_CTLA4_CD80*x_syn_T_C1_CTLA4)*x_syn_T_C1_CD80_CTLA4)-((2*koff_CTLA4_CD80)*x_syn_T_C1_CTLA4_CD80_CTLA4))   # R123 syn_T_C1.CTLA4 + syn_T_C1.CD80_CTLA4 <-> syn_T_C1.CTLA4_CD80_CTLA4: binding and unbinding of CTLA4 and CD80-CTLA4 in synapse
    v124 <- (((kon_CTLA4_CD80*x_syn_T_C1_CD80)*x_syn_T_C1_CTLA4_CD80_CTLA4)-(koff_CTLA4_CD80*x_syn_T_C1_CD80_CTLA4_CD80_CTLA4))   # R124 syn_T_C1.CD80 + syn_T_C1.CTLA4_CD80_CTLA4 <-> syn_T_C1.CD80_CTLA4_CD80_CTLA4: binding and unbinding of CD80 and CTLA4-CD80-CTLA4 in synapse
    v125 <- (((kon_CTLA4_CD80*x_syn_T_C1_CD80_CTLA4)*x_syn_T_C1_CD80)-((2*koff_CTLA4_CD80)*x_syn_T_C1_CD80_CTLA4_CD80))   # R125 syn_T_C1.CD80_CTLA4 + syn_T_C1.CD80 <-> syn_T_C1.CD80_CTLA4_CD80: binding and unbinding of CD80-CTLA4 and CD80 in synapse
    v126 <- (((kon_CTLA4_CD80*x_syn_T_C1_CTLA4)*x_syn_T_C1_CD80_CTLA4_CD80)-(koff_CTLA4_CD80*x_syn_T_C1_CD80_CTLA4_CD80_CTLA4))   # R126 syn_T_C1.CTLA4 + syn_T_C1.CD80_CTLA4_CD80 <-> syn_T_C1.CD80_CTLA4_CD80_CTLA4: binding and unbinding of CTLA4 and CD80-CTLA4-CD80 in synapse
    v127 <- ((((2*kon_CTLA4_CD86)*x_syn_T_C1_CTLA4)*x_syn_T_C1_CD86)-(koff_CTLA4_CD86*x_syn_T_C1_CD86_CTLA4))   # R127 syn_T_C1.CTLA4 + syn_T_C1.CD86 <-> syn_T_C1.CD86_CTLA4: binding and unbinding of CTLA4 and CD86 in synapse
    v128 <- (((kon_CTLA4_CD86*x_syn_T_C1_CD86_CTLA4)*x_syn_T_C1_CD86)-((2*koff_CTLA4_CD86)*x_syn_T_C1_CD86_CTLA4_CD86))   # R128 syn_T_C1.CD86_CTLA4 + syn_T_C1.CD86 <-> syn_T_C1.CD86_CTLA4_CD86: binding and unbinding of CD86-CTLA4 and CD86 in synapse
    v129 <- (((4*kon_CTLA4_ipi)*((x_syn_T_C1_CTLA4*x_V_T_ipi)/gamma_T_ipi))-(koff_CTLA4_ipi*x_syn_T_C1_CTLA4_ipi))   # R129 syn_T_C1.CTLA4 <-> syn_T_C1.CTLA4_ipi: binding and unbinding of CTLA4 to Ipi on T surface in synapse
    v130 <- (((Chi_CTLA4_ipi*kon_CTLA4_ipi)*(x_syn_T_C1_CTLA4*x_syn_T_C1_CTLA4_ipi))-((2*koff_CTLA4_ipi)*x_syn_T_C1_CTLA4_ipi_CTLA4))   # R130 syn_T_C1.CTLA4_ipi + syn_T_C1.CTLA4 <-> syn_T_C1.CTLA4_ipi_CTLA4: binding and unbinding of CTLA4 to Ipi on T surface in synapse
    v131 <- ((((2*kon_CD80_PDL1)*x_syn_T_C1_CD80)*x_syn_T_C1_TPDL1)-(koff_CD80_PDL1*x_syn_T_C1_TPDL1_CD80))   # R131 syn_T_C1.CD80 + syn_T_C1.TPDL1 <-> syn_T_C1.TPDL1_CD80: binding and unbinding of CD80 and PDL1 in synapse
    v132 <- (((kon_CD80_PDL1*x_syn_T_C1_TPDL1_CD80)*x_syn_T_C1_TPDL1)-((2*koff_CD80_PDL1)*x_syn_T_C1_TPDL1_CD80_TPDL1))   # R132 syn_T_C1.TPDL1_CD80 + syn_T_C1.TPDL1 <-> syn_T_C1.TPDL1_CD80_TPDL1: binding and unbinding of PDL1-CD80 and PDL1 in synapse
    v133 <- (((2*kon_PDL1_atezo)*((x_syn_T_C1_TPDL1*x_V_T_atezo)/gamma_T_atezo))-(koff_PDL1_atezo*x_syn_T_C1_TPDL1_atezo))   # R133 syn_T_C1.TPDL1 <-> syn_T_C1.TPDL1_atezo: binding and unbinding of PDL1 to atezo on C1 surface in synapse
    v134 <- (((Chi_PDL1_atezo*kon_PDL1_atezo)*(x_syn_T_C1_TPDL1*x_syn_T_C1_TPDL1_atezo))-((2*koff_PDL1_atezo)*x_syn_T_C1_TPDL1_atezo_TPDL1))   # R134 syn_T_C1.TPDL1_atezo + syn_T_C1.TPDL1 <-> syn_T_C1.TPDL1_atezo_TPDL1: binding and unbinding of PDL1 to PDL1-atezo on C1 surface in synapse
    v135 <- (((k_out_PDL1*x_V_T_IFNg)/(x_V_T_IFNg+a_IFNg_50_ind))*(1-(x_syn_T_APC_PDL1_total/((APC_PDL1_base*r_PDL1_IFNg)/A_cell))))   # R135 null -> syn_T_APC.PDL1: Translocation of PDL1 between cell surface and cytoplasm
    v136 <- ((((k_out_PDL1*r_PDL2APC)*x_V_T_IFNg)/(x_V_T_IFNg+a_IFNg_50_ind))*(1-(x_syn_T_APC_PDL2_total/(((APC_PDL1_base*r_PDL1_IFNg)/A_cell)*r_PDL2APC))))   # R136 null -> syn_T_APC.PDL2: Translocation of PDL2 between cell surface and cytoplasm
    v137 <- (k_in_PDL1*((APC_PDL1_base/A_cell)-x_syn_T_APC_PDL1_total))   # R137 null -> syn_T_APC.PDL1: Translocation of PDL1 between cell surface and cytoplasm
    v138 <- (k_in_PDL1*(((APC_PDL1_base/A_cell)*r_PDL2APC)-x_syn_T_APC_PDL2_total))   # R138 null -> syn_T_APC.PDL2: Translocation of PDL2 between cell surface and cytoplasm
    v139 <- (((kon_PD1_PDL1*x_syn_T_APC_PD1)*x_syn_T_APC_PDL1)-(koff_PD1_PDL1*x_syn_T_APC_PD1_PDL1))   # R139 syn_T_APC.PD1 + syn_T_APC.PDL1 <-> syn_T_APC.PD1_PDL1: binding and unbinding of PD1 PDL1 in synapse
    v140 <- (((kon_PD1_PDL2*x_syn_T_APC_PD1)*x_syn_T_APC_PDL2)-(koff_PD1_PDL2*x_syn_T_APC_PD1_PDL2))   # R140 syn_T_APC.PD1 + syn_T_APC.PDL2 <-> syn_T_APC.PD1_PDL2: binding and unbinding of PD1 PDL2 in synapse
    v141 <- (((2*kon_PD1_nivo)*((x_syn_T_APC_PD1*x_V_LN_nivo)/gamma_LN_nivo))-(koff_PD1_nivo*x_syn_T_APC_PD1_nivo))   # R141 syn_T_APC.PD1 <-> syn_T_APC.PD1_nivo: binding and unbinding of PD1 to Nivo on T surface in synapse
    v142 <- (((Chi_PD1_nivo*kon_PD1_nivo)*(x_syn_T_APC_PD1*x_syn_T_APC_PD1_nivo))-((2*koff_PD1_nivo)*x_syn_T_APC_PD1_nivo_PD1))   # R142 syn_T_APC.PD1_nivo + syn_T_APC.PD1 <-> syn_T_APC.PD1_nivo_PD1: binding and unbinding of PD1 to PD1-Nivo on T surface in synapse
    v143 <- (((2*kon_PDL1_atezo)*((x_syn_T_APC_PDL1*x_V_LN_atezo)/gamma_LN_atezo))-(koff_PDL1_atezo*x_syn_T_APC_PDL1_atezo))   # R143 syn_T_APC.PDL1 <-> syn_T_APC.PDL1_atezo: binding and unbinding of PDL1 to atezo on APC surface in synapse
    v144 <- (((Chi_PDL1_atezo*kon_PDL1_atezo)*(x_syn_T_APC_PDL1*x_syn_T_APC_PDL1_atezo))-((2*koff_PDL1_atezo)*x_syn_T_APC_PDL1_atezo_PDL1))   # R144 syn_T_APC.PDL1_atezo + syn_T_APC.PDL1 <-> syn_T_APC.PDL1_atezo_PDL1: binding and unbinding of PDL1 to PDL1-atezo on APC surface in synapse
    v145 <- ((((2*kon_CD28_CD80)*x_syn_T_APC_CD28)*x_syn_T_APC_CD80)-(koff_CD28_CD80*x_syn_T_APC_CD28_CD80))   # R145 syn_T_APC.CD28 + syn_T_APC.CD80 <-> syn_T_APC.CD28_CD80: binding and unbinding of CD28 and CD80 in synapse
    v146 <- (((kon_CD28_CD80*x_syn_T_APC_CD28)*x_syn_T_APC_CD28_CD80)-((2*koff_CD28_CD80)*x_syn_T_APC_CD28_CD80_CD28))   # R146 syn_T_APC.CD28_CD80 + syn_T_APC.CD28 <-> syn_T_APC.CD28_CD80_CD28: binding and unbinding of CD28-CD80 and CD28 in synapse
    v147 <- (((kon_CD28_CD86*x_syn_T_APC_CD28)*x_syn_T_APC_CD86)-(koff_CD28_CD86*x_syn_T_APC_CD28_CD86))   # R147 syn_T_APC.CD28 + syn_T_APC.CD86 <-> syn_T_APC.CD28_CD86: binding and unbinding of CD28 and CD86 in synapse
    v148 <- ((((4*kon_CTLA4_CD80)*x_syn_T_APC_CTLA4)*x_syn_T_APC_CD80)-(koff_CTLA4_CD80*x_syn_T_APC_CD80_CTLA4))   # R148 syn_T_APC.CTLA4 + syn_T_APC.CD80 <-> syn_T_APC.CD80_CTLA4: binding and unbinding of CTLA4 and CD80 in synapse
    v149 <- (((kon_CTLA4_CD80*x_syn_T_APC_CTLA4)*x_syn_T_APC_CD80_CTLA4)-((2*koff_CTLA4_CD80)*x_syn_T_APC_CTLA4_CD80_CTLA4))   # R149 syn_T_APC.CTLA4 + syn_T_APC.CD80_CTLA4 <-> syn_T_APC.CTLA4_CD80_CTLA4: binding and unbinding of CTLA4 and CD80-CTLA4 in synapse
    v150 <- (((kon_CTLA4_CD80*x_syn_T_APC_CD80)*x_syn_T_APC_CTLA4_CD80_CTLA4)-(koff_CTLA4_CD80*x_syn_T_APC_CD80_CTLA4_CD80_CTLA4))   # R150 syn_T_APC.CD80 + syn_T_APC.CTLA4_CD80_CTLA4 <-> syn_T_APC.CD80_CTLA4_CD80_CTLA4: binding and unbinding of CD80 and CTLA4-CD80-CTLA4 in synapse
    v151 <- (((kon_CTLA4_CD80*x_syn_T_APC_CD80_CTLA4)*x_syn_T_APC_CD80)-((2*koff_CTLA4_CD80)*x_syn_T_APC_CD80_CTLA4_CD80))   # R151 syn_T_APC.CD80_CTLA4 + syn_T_APC.CD80 <-> syn_T_APC.CD80_CTLA4_CD80: binding and unbinding of CD80-CTLA4 and CD80 in synapse
    v152 <- (((kon_CTLA4_CD80*x_syn_T_APC_CTLA4)*x_syn_T_APC_CD80_CTLA4_CD80)-(koff_CTLA4_CD80*x_syn_T_APC_CD80_CTLA4_CD80_CTLA4))   # R152 syn_T_APC.CTLA4 + syn_T_APC.CD80_CTLA4_CD80 <-> syn_T_APC.CD80_CTLA4_CD80_CTLA4: binding and unbinding of CTLA4 and CD80-CTLA4-CD80 in synapse
    v153 <- ((((2*kon_CTLA4_CD86)*x_syn_T_APC_CTLA4)*x_syn_T_APC_CD86)-(koff_CTLA4_CD86*x_syn_T_APC_CD86_CTLA4))   # R153 syn_T_APC.CTLA4 + syn_T_APC.CD86 <-> syn_T_APC.CD86_CTLA4: binding and unbinding of CTLA4 and CD86 in synapse
    v154 <- (((kon_CTLA4_CD86*x_syn_T_APC_CD86_CTLA4)*x_syn_T_APC_CD86)-((2*koff_CTLA4_CD86)*x_syn_T_APC_CD86_CTLA4_CD86))   # R154 syn_T_APC.CD86_CTLA4 + syn_T_APC.CD86 <-> syn_T_APC.CD86_CTLA4_CD86: binding and unbinding of CD86-CTLA4 and CD86 in synapse
    v155 <- (((4*kon_CTLA4_ipi)*((x_syn_T_APC_CTLA4*x_V_LN_ipi)/gamma_LN_ipi))-(koff_CTLA4_ipi*x_syn_T_APC_CTLA4_ipi))   # R155 syn_T_APC.CTLA4 <-> syn_T_APC.CTLA4_ipi: binding and unbinding of CTLA4 to Ipi on T surface in synapse
    v156 <- (((Chi_CTLA4_ipi*kon_CTLA4_ipi)*(x_syn_T_APC_CTLA4*x_syn_T_APC_CTLA4_ipi))-((2*koff_CTLA4_ipi)*x_syn_T_APC_CTLA4_ipi_CTLA4))   # R156 syn_T_APC.CTLA4_ipi + syn_T_APC.CTLA4 <-> syn_T_APC.CTLA4_ipi_CTLA4: binding and unbinding of CTLA4 to Ipi on T surface in synapse
    v157 <- ((((2*kon_CD80_PDL1)*x_syn_T_APC_CD80)*x_syn_T_APC_TPDL1)-(koff_CD80_PDL1*x_syn_T_APC_TPDL1_CD80))   # R157 syn_T_APC.CD80 + syn_T_APC.TPDL1 <-> syn_T_APC.TPDL1_CD80: binding and unbinding of CD80 and PDL1 in synapse
    v158 <- (((kon_CD80_PDL1*x_syn_T_APC_TPDL1_CD80)*x_syn_T_APC_TPDL1)-((2*koff_CD80_PDL1)*x_syn_T_APC_TPDL1_CD80_TPDL1))   # R158 syn_T_APC.TPDL1_CD80 + syn_T_APC.TPDL1 <-> syn_T_APC.TPDL1_CD80_TPDL1: binding and unbinding of PDL1-CD80 and PDL1 in synapse
    v159 <- (((2*kon_PDL1_atezo)*((x_syn_T_APC_TPDL1*x_V_LN_atezo)/gamma_LN_atezo))-(koff_PDL1_atezo*x_syn_T_APC_TPDL1_atezo))   # R159 syn_T_APC.TPDL1 <-> syn_T_APC.TPDL1_atezo: binding and unbinding of PDL1 to atezo on APC surface in synapse
    v160 <- (((Chi_PDL1_atezo*kon_PDL1_atezo)*(x_syn_T_APC_TPDL1*x_syn_T_APC_TPDL1_atezo))-((2*koff_PDL1_atezo)*x_syn_T_APC_TPDL1_atezo_TPDL1))   # R160 syn_T_APC.TPDL1_atezo + syn_T_APC.TPDL1 <-> syn_T_APC.TPDL1_atezo_TPDL1: binding and unbinding of PDL1 to PDL1-atezo on APC surface in synapse
    v161 <- ((kon_CTLA4_ipi*((x_V_T_Treg_CTLA4*x_V_T_ipi)/gamma_C_ipi))-(koff_CTLA4_ipi*x_V_T_Treg_CTLA4_ipi))   # R161 V_T.Treg_CTLA4 <-> V_T.Treg_CTLA4_ipi: binding and unbinding of CTLA4 to Ipi on Treg surface in Tumor
    v162 <- ((((Chi_CTLA4_ipi*kon_CTLA4_ipi)*(x_V_T_Treg_CTLA4*x_V_T_Treg_CTLA4_ipi))/A_Tcell)-(koff_CTLA4_ipi*x_V_T_Treg_CTLA4_ipi_CTLA4))   # R162 V_T.Treg_CTLA4_ipi + V_T.Treg_CTLA4 <-> V_T.Treg_CTLA4_ipi_CTLA4: binding and unbinding of CTLA4 to CTLA4-Ipi on Treg surface in Tumor
    v163 <- ((kon_CTLA4_ipi*((x_V_P_Treg_CTLA4*x_V_P_ipi)/gamma_P_ipi))-(koff_CTLA4_ipi*x_V_P_Treg_CTLA4_ipi))   # R163 V_P.Treg_CTLA4 <-> V_P.Treg_CTLA4_ipi: binding and unbinding of CTLA4 to Ipi on Treg surface in Peripheral compartment
    v164 <- ((((Chi_CTLA4_ipi*kon_CTLA4_ipi)*(x_V_P_Treg_CTLA4*x_V_P_Treg_CTLA4_ipi))/A_Tcell)-(koff_CTLA4_ipi*x_V_P_Treg_CTLA4_ipi_CTLA4))   # R164 V_P.Treg_CTLA4_ipi + V_P.Treg_CTLA4 <-> V_P.Treg_CTLA4_ipi_CTLA4: binding and unbinding of CTLA4 to CTLA4-Ipi on Treg surface in Peripheral compartment
    v165 <- ((k_CTLA4_ADCC*x_V_T_T0)*H_Treg_T)   # R165 V_T.T0 -> null: CTLA4 ADCC in the tumor compartment
    v166 <- ((k_CTLA4_ADCC*x_V_P_T0)*H_Treg_P)   # R166 V_P.T0 -> null: CTLA4 ADCC in the peripheral compartment
    v167 <- ((k_rec_MDSC*((MDSC_max*vol_V_T)-x_V_T_MDSC))*(x_V_T_CCL2/(a_EC50_CCL2_rec+x_V_T_CCL2)))   # R167 null -> V_T.MDSC: MDSC recruitment into the tumor compartment
    v168 <- (k_brec_MDSC*((MDSC_max*vol_V_T)-x_V_T_MDSC))   # R168 null -> V_T.MDSC: Baseline MDSC Migration into the tumor compartment
    v169 <- (kd_MDSC*x_V_T_MDSC)   # R169 V_T.MDSC -> null: MDSC death in the tumor compartment
    v170 <- (k_deg_CCL2*x_V_T_CCL2)   # R170 V_T.CCL2 -> null: CCL2 degradation
    v171 <- (k_deg_NO*x_V_T_NO)   # R171 V_T.NO -> null: NO degradation
    v172 <- (k_deg_ArgI*x_V_T_ArgI)   # R172 V_T.ArgI -> null: ArgI degradation
    v173 <- ((k_sec_CCL2*C_total)*(1-(x_V_T_ENT/(x_V_T_ENT+a_IC50_ENT_CCL2))))   # R173 null -> V_T.CCL2: CCL2 secretion by MDSCs and cancer cells
    v174 <- ((k_sec_NO*x_V_T_MDSC)*(1-(x_V_T_ENT/(x_V_T_ENT+a_IC50_ENT_NO))))   # R174 null -> V_T.NO: NO secretion from MDSC
    v175 <- ((k_sec_ArgI*x_V_T_MDSC)*(1-(x_V_T_ENT/(x_V_T_ENT+a_IC50_ENT_ArgI))))   # R175 null -> V_T.ArgI: ArgI secretion from MDSC
    v176 <- (k_dose2*x_V_C_Dose2)   # R176 V_C.Dose2 -> V_C.ENT_GI: first-order absorption of ENT into buccal compartment
    v177 <- (k_a1_ENT*x_V_C_ENT_Buccal)   # R177 V_C.ENT_Buccal -> V_C.ENT: Buccal absorption of ENTinto central comparment
    v178 <- (k_a2_ENT*x_V_C_ENT_GI)   # R178 V_C.ENT_GI -> V_C.ENT: Gastro-intestinal absorption of ENTinto central comparment
    v179 <- ((k_cln_ENT*x_V_C_ENT)/(x_V_C_ENT+a_Kc_ENT))   # R179 V_C.ENT -> null: ENT Non-linear clearance from central compartment
    v180 <- (q_P_ENT*((x_V_C_ENT/gamma_C_ENT)-(x_V_P_ENT/gamma_P_ENT)))   # R180 V_C.ENT <-> V_P.ENT: ENT diffusive transport to peripheral compartment
    v181 <- (q_T_ENT*((x_V_C_ENT/gamma_C_ENT)-(x_V_T_ENT/gamma_T_ENT)))   # R181 V_C.ENT <-> V_T.ENT: ENT diffusive transport to tumor compartment
    v182 <- (q_LN_ENT*((x_V_C_ENT/gamma_C_ENT)-(x_V_LN_ENT/gamma_LN_ENT)))   # R182 V_C.ENT <-> V_LN.ENT: ENT diffusive transport to lymph node compartment
    v183 <- ((q_LD_ENT*x_V_T_ENT)/gamma_T_ENT)   # R183 V_T.ENT -> V_LN.ENT: ENT convective transport from tumor to lymph node
    v184 <- ((q_LD_ENT*x_V_LN_ENT)/gamma_LN_ENT)   # R184 V_LN.ENT -> V_C.ENT: ENT convective transport from lymph node to central
    v185 <- (k_cl_ENT*x_V_C_ENT)   # R185 V_C.ENT -> null: Drug clearance from central compartment
    v186 <- (((Vmt/(x_V_1_NabP+a_Kt))*x_V_1_NabP)-((Vmt/(x_V_2_NabP+a_Kt))*x_V_2_NabP))   # R186 V_1.NabP <-> V_2.NabP: Intercompartmental distribution of nab-paclitaxel between V1 and V2 compartment
    v187 <- ((Q2*x_V_1_NabP)-(Q2*x_V_3_NabP))   # R187 V_1.NabP <-> V_3.NabP: Intercompartmental clearance of nab-paclitaxel between V1 and V3 compartment
    v188 <- ((Vmcl/(x_V_1_NabP+a_Kcl))*x_V_1_NabP)   # R188 V_1.NabP -> null: Clearance of nab-paclitaxel from V1 compartment
    v189 <- ((((k_C_nabp*x_V_T_C1)*(x_V_T_NabP/(x_V_T_NabP+a_IC50_nabp)))*min(C_total, Kc_cyt))/C_total)   # R189 V_T.C1 -> V_T.C_x: Cancer cell death by nab-paclitaxel
    v190 <- ((k_C_resist*x_V_T_C1)*H_TGF_CTL)   # R190 V_T.C1 -> V_T.C2: Cancer cell resistance to nab-paclitaxel
    v191 <- ((((k_C_nabp*x_V_T_C2)*(x_V_T_NabP/(x_V_T_NabP+(a_IC50_nabp*r_resist))))*min(C_total, Kc_cyt))/C_total)   # R191 V_T.C2 -> V_T.C_x: Resistant cancer cell death by nab-paclitaxel
    v192 <- (k_c_vas*C_total)   # R192 null -> V_T.c_vas: Baseline secretion/degradation of angiogenic factors
    v193 <- (((k_c_nabp*C_total)*x_V_T_NabP)/(x_V_T_NabP+a_IC50_nabp_vas))   # R193 null -> V_T.c_vas: Angiogenic factor release in response to nab-paclitaxel
    v194 <- (k_c_deg*x_V_T_c_vas)   # R194 V_T.c_vas -> null: Baseline secretion/degradation of angiogenic factors
    v195 <- (((k_K_g*C_total)*x_V_T_c_vas)/(x_V_T_c_vas+a_IC50_vas))   # R195 null -> V_T.K: Resistant cancer cell death by nab-paclitaxel
    v196 <- ((k_K_d*x_V_T_K)*((((C_total/a_cell)*2.57e-06)^(1.0/3.0))^2))   # R196 V_T.K -> null: Resistant cancer cell death by nab-paclitaxel
    v197 <- ((k_endo_nabp*x_V_T_K)*x_V_T_NabP)   # R197 V_T.K -> null: Endothelial cell death by nab-paclitaxel
    v198 <- (((k_Th_act*H_APCh)*H_P1)*x_V_LN_nT0)   # R198 V_LN.nT0 -> null: Naive T cell activation
    v199 <- ((((k_Th_act*H_APCh)*H_P1)*x_V_LN_nT0)*n_T1_clones)   # R199 null -> V_LN.aTh: Naive T cell activation
    v200 <- ((k_T0_pro/N_aTh)*x_V_LN_aTh)   # R200 V_LN.aTh -> null: aTh cell proliferation
    v201 <- (((k_T0_pro/N_aTh)*(2^N_aTh))*x_V_LN_aTh)   # R201 null -> V_LN.Th: aTh cell proliferation
    v202 <- ((((k_reg*x_V_T_Th)*x_V_T_ArgI)/(a_EC50_ArgI_Treg+x_V_T_ArgI))*H_TGF_reg)   # R202 V_T.Th -> V_T.T0: Differentiation between Treg and Th Cells
    v203 <- (k_T0_death*x_V_C_Th)   # R203 V_C.Th -> null: T cell death in the central compartment
    v204 <- (k_T0_death*x_V_P_Th)   # R204 V_P.Th -> null: T cell death in the peripheral compartment
    v205 <- (k_T0_death*x_V_LN_Th)   # R205 V_LN.Th -> null: T cell death in the lymph node compartment
    v206 <- (k_T0_death*x_V_T_Th)   # R206 V_T.Th -> V_T.Th_exh: T cell death in the tumor compartment
    v207 <- (q_T0_P_in*x_V_C_Th)   # R207 V_C.Th -> V_P.Th: T cell transport into the peripheral compartment
    v208 <- (q_T0_P_out*x_V_P_Th)   # R208 V_P.Th -> V_C.Th: T cell transport out of the peripheral compartment
    v209 <- (((q_T0_T_in*vol_V_T)*x_V_C_Th)*((C_total^2)/((C_total^2)+a_Kc_rec)))   # R209 V_C.Th -> V_T.Th: T cell transport into the tumor compartment
    v210 <- (q_T0_LN_out*x_V_LN_Th)   # R210 V_LN.Th -> V_C.Th: T cell transport out of the lymph node compartment
    v211 <- (k_IL2_sec*x_V_LN_aTh)   # R211 null -> V_LN.IL2: IL2 secretion from activated T cells
    v212 <- ((k_TGF_deg*(a_TGFbase-x_V_T_TGF))*vol_V_T)   # R212 null -> V_T.TGF: TGFb secretion by triple-negative breast cancer cells
    v213 <- (k_TGF_Tsec*x_V_T_T0)   # R213 null -> V_T.TGF: TGF secretion from activated T cells in tumor
    v214 <- (k_IFNg_sec*x_V_T_Th)   # R214 null -> V_T.IFNg: IFNg secretion from T helper cells in tumor
    v215 <- (k_IFNg_deg*x_V_T_IFNg)   # R215 V_T.IFNg -> null: IFNg degradtion in tumor
    v216 <- ((((k_C_carb*x_V_T_C1)*(x_V_T_carb/(x_V_T_carb+a_IC50_carb)))*min(C_total, Kc_carb))/C_total)   # R216 V_T.C1 -> V_T.C_x: Cancer cell death by nab-paclitaxel

    # ---------- ODEs (states are AMOUNTS in the canonical units) ----------
    d/dt(q_V_C_nT0) <- v8-v12-v15+v16-v17+v18
    d/dt(q_V_C_T0) <- -v23-v27+v28-v29+v30
    d/dt(q_V_C_nT1) <- v34-v38-v41+v42-v43+v44
    d/dt(q_V_C_T1) <- -v49-v55+v56-v57+v58
    d/dt(q_V_C_nivo) <- -v91-v92-v93+v95*vol_V_T-v96
    d/dt(q_V_C_atezo) <- -v97-v98-v99+v101*vol_V_T-v102
    d/dt(q_V_C_ipi) <- -v103-v104-v105+v107*vol_V_T-v108
    d/dt(q_V_C_ENT) <- v177*m_V_C+v178*m_V_C-v179*m_V_C-v180-v181-v182+v184*vol_V_T-v185*m_V_C
    d/dt(q_V_C_ENT_Buccal) <- -v177*m_V_C
    d/dt(q_V_C_ENT_GI) <- v176*m_V_C-v178*m_V_C
    d/dt(q_V_C_Dose2) <- -v176*m_V_C
    d/dt(q_V_C_Th) <- -v203-v207+v208-v209+v210
    d/dt(q_V_P_nT0) <- v9-v11+v15-v16
    d/dt(q_V_P_T0) <- -v24+v27-v28-v166
    d/dt(q_V_P_nT1) <- v35-v37+v41-v42
    d/dt(q_V_P_T1) <- -v50+v55-v56
    d/dt(q_V_P_nivo) <- v91
    d/dt(q_V_P_atezo) <- v97
    d/dt(q_V_P_ipi) <- v103
    d/dt(q_V_P_Treg_CTLA4) <- -v163*NAVG_nmol-v164*NAVG_nmol
    d/dt(q_V_P_Treg_CTLA4_ipi) <- v163*NAVG_nmol-v164*NAVG_nmol
    d/dt(q_V_P_Treg_CTLA4_ipi_CTLA4) <- v164*NAVG_nmol
    d/dt(q_V_P_ENT) <- v180
    d/dt(q_V_P_Th) <- -v204+v207-v208
    d/dt(q_V_T_C_x) <- -v1+v5+v7+v60+v61+v189+v191+v216
    d/dt(q_V_T_T1_exh) <- -v2+v52+v53+v54
    d/dt(q_V_T_Th_exh) <- -v3+v206
    dq_V_T_C1 <- v4-v5-v60-v189-v190-v216
    if (q_V_T_C1 < 0.5 && dq_V_T_C1 < 0) {
      dq_V_T_C1 <- 0   # Table S7 event 1: floor at 0.01 cell
    }
    d/dt(q_V_T_C1) <- dq_V_T_C1
    dq_V_T_C2 <- v6-v7-v61+v190-v191
    if (q_V_T_C2 < 0.5 && dq_V_T_C2 < 0) {
      dq_V_T_C2 <- 0   # Table S7 event 2: floor at 0.01 cell
    }
    d/dt(q_V_T_C2) <- dq_V_T_C2
    d/dt(q_V_T_nT0) <- -v14
    d/dt(q_V_T_T0) <- -v26+v29-v165+v202
    d/dt(q_V_T_nT1) <- -v40
    d/dt(q_V_T_T1) <- -v52-v53-v54+v57
    d/dt(q_V_T_APC) <- v62-v64
    d/dt(q_V_T_mAPC) <- v64-v65-v66
    d/dt(q_V_T_c) <- v68*vol_V_T+v69
    d/dt(q_V_T_P0) <- v71-v72-v73
    d/dt(q_V_T_P1) <- v81-v82-v83
    d/dt(q_V_T_nivo) <- v92-v94*vol_V_T
    d/dt(q_V_T_atezo) <- v98-v100*vol_V_T
    d/dt(q_V_T_ipi) <- v104-v106*vol_V_T
    d/dt(q_V_T_Treg_CTLA4) <- -v161*NAVG_nmol-v162*NAVG_nmol
    d/dt(q_V_T_Treg_CTLA4_ipi) <- v161*NAVG_nmol-v162*NAVG_nmol
    d/dt(q_V_T_Treg_CTLA4_ipi_CTLA4) <- v162*NAVG_nmol
    d/dt(q_V_T_MDSC) <- v167+v168-v169
    d/dt(q_V_T_CCL2) <- -v170*vol_V_T+v173
    d/dt(q_V_T_NO) <- -v171*vol_V_T+v174
    d/dt(q_V_T_ArgI) <- -v172+v175/vol_V_T
    d/dt(q_V_T_ENT) <- v181-v183*vol_V_T
    dq_V_T_K <- v195-v196-v197
    if (q_V_T_K < 0.5 && dq_V_T_K < 0) {
      dq_V_T_K <- 0   # Table S7 event 3: floor at 0.01 cell
    }
    d/dt(q_V_T_K) <- dq_V_T_K
    d/dt(q_V_T_c_vas) <- v192+v193-v194*vol_V_T
    d/dt(q_V_T_Th) <- -v202-v206+v209
    d/dt(q_V_T_TGF) <- v212+v213
    d/dt(q_V_T_IFNg) <- v214-v215*vol_V_T
    d/dt(q_V_LN_nT0) <- v10-v13+v17-v18-v19-v198
    d/dt(q_V_LN_aT0) <- v20-v21
    d/dt(q_V_LN_T0) <- v22-v25-v30
    d/dt(q_V_LN_IL2) <- -v31-v32-v33+v59+v211
    d/dt(q_V_LN_nT1) <- v36-v39+v43-v44-v45
    d/dt(q_V_LN_aT1) <- v46-v47
    d/dt(q_V_LN_T1) <- v48-v51-v58
    d/dt(q_V_LN_APC) <- v63
    d/dt(q_V_LN_mAPC) <- v65-v67
    d/dt(q_V_LN_nivo) <- v93+v94*vol_V_T-v95*vol_V_T
    d/dt(q_V_LN_atezo) <- v99+v100*vol_V_T-v101*vol_V_T
    d/dt(q_V_LN_ipi) <- v105+v106*vol_V_T-v107*vol_V_T
    d/dt(q_V_LN_ENT) <- v182+v183*vol_V_T-v184*vol_V_T
    d/dt(q_V_LN_aTh) <- v199-v200
    d/dt(q_V_LN_Th) <- v201-v205-v210
    d/dt(q_V_e_P0) <- v74-v75
    d/dt(q_V_e_p0) <- v75-v76-v77+v78
    d/dt(q_V_e_P1) <- v84-v85
    d/dt(q_V_e_p1) <- v85-v86-v87+v88
    d/dt(q_A_e_M1) <- -v70*NAVG_nmol-v77*NAVG_nmol+v78*NAVG_nmol-v87*NAVG_nmol+v88*NAVG_nmol
    d/dt(q_A_e_M1p0) <- v77*NAVG_nmol-v78*NAVG_nmol-v80*NAVG_nmol
    d/dt(q_A_e_M1p1) <- v87*NAVG_nmol-v88*NAVG_nmol-v90*NAVG_nmol
    d/dt(q_A_s_M1) <- v70*NAVG_nmol+v79*NAVG_nmol+v89*NAVG_nmol
    d/dt(q_A_s_M1p0) <- -v79*NAVG_nmol+v80*NAVG_nmol
    d/dt(q_A_s_M1p1) <- -v89*NAVG_nmol+v90*NAVG_nmol
    d/dt(q_syn_T_C1_PD1_PDL1) <- v113*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_PD1_PDL2) <- v114*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_PD1) <- -v113*m_syn_T_C1*NAVG_nmol-v114*m_syn_T_C1*NAVG_nmol-v115*m_syn_T_C1*NAVG_nmol-v116*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_PDL1) <- v109*NAVG_nmol+v111*m_syn_T_C1*NAVG_nmol-v113*m_syn_T_C1*NAVG_nmol-v117*m_syn_T_C1*NAVG_nmol-v118*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_PDL2) <- v110*NAVG_nmol+v112*m_syn_T_C1*NAVG_nmol-v114*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_PD1_nivo) <- v115*m_syn_T_C1*NAVG_nmol-v116*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_PD1_nivo_PD1) <- v116*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_PDL1_atezo) <- v117*m_syn_T_C1*NAVG_nmol-v118*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_PDL1_atezo_PDL1) <- v118*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_TPDL1) <- -v131*m_syn_T_C1*NAVG_nmol-v132*m_syn_T_C1*NAVG_nmol-v133*m_syn_T_C1*NAVG_nmol-v134*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_TPDL1_atezo) <- v133*m_syn_T_C1*NAVG_nmol-v134*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_TPDL1_atezo_TPDL1) <- v134*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_CD28_CD80) <- v119*m_syn_T_C1*NAVG_nmol-v120*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_CD28_CD80_CD28) <- v120*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_CD28_CD86) <- v121*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_CD80_CTLA4) <- v122*m_syn_T_C1*NAVG_nmol-v123*m_syn_T_C1*NAVG_nmol-v125*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_CD80_CTLA4_CD80) <- v125*m_syn_T_C1*NAVG_nmol-v126*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_CTLA4_CD80_CTLA4) <- v123*m_syn_T_C1*NAVG_nmol-v124*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_CD80_CTLA4_CD80_CTLA4) <- v124*m_syn_T_C1*NAVG_nmol+v126*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_CD86_CTLA4) <- v127*m_syn_T_C1*NAVG_nmol-v128*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_CD86_CTLA4_CD86) <- v128*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_TPDL1_CD80) <- v131*m_syn_T_C1*NAVG_nmol-v132*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_TPDL1_CD80_TPDL1) <- v132*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_CD28) <- -v119*m_syn_T_C1*NAVG_nmol-v120*m_syn_T_C1*NAVG_nmol-v121*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_CTLA4) <- -v122*m_syn_T_C1*NAVG_nmol-v123*m_syn_T_C1*NAVG_nmol-v126*m_syn_T_C1*NAVG_nmol-v127*m_syn_T_C1*NAVG_nmol-v129*m_syn_T_C1*NAVG_nmol-v130*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_CD80) <- -v119*m_syn_T_C1*NAVG_nmol-v122*m_syn_T_C1*NAVG_nmol-v124*m_syn_T_C1*NAVG_nmol-v125*m_syn_T_C1*NAVG_nmol-v131*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_CD86) <- -v121*m_syn_T_C1*NAVG_nmol-v127*m_syn_T_C1*NAVG_nmol-v128*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_CTLA4_ipi) <- v129*m_syn_T_C1*NAVG_nmol-v130*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_C1_CTLA4_ipi_CTLA4) <- v130*m_syn_T_C1*NAVG_nmol
    d/dt(q_syn_T_APC_PD1_PDL1) <- v139*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_PD1_PDL2) <- v140*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_PD1) <- -v139*m_syn_T_APC*NAVG_nmol-v140*m_syn_T_APC*NAVG_nmol-v141*m_syn_T_APC*NAVG_nmol-v142*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_PDL1) <- v135*NAVG_nmol+v137*m_syn_T_APC*NAVG_nmol-v139*m_syn_T_APC*NAVG_nmol-v143*m_syn_T_APC*NAVG_nmol-v144*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_PDL2) <- v136*NAVG_nmol+v138*m_syn_T_APC*NAVG_nmol-v140*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_PD1_nivo) <- v141*m_syn_T_APC*NAVG_nmol-v142*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_PD1_nivo_PD1) <- v142*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_PDL1_atezo) <- v143*m_syn_T_APC*NAVG_nmol-v144*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_PDL1_atezo_PDL1) <- v144*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_TPDL1) <- -v157*m_syn_T_APC*NAVG_nmol-v158*m_syn_T_APC*NAVG_nmol-v159*m_syn_T_APC*NAVG_nmol-v160*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_TPDL1_atezo) <- v159*m_syn_T_APC*NAVG_nmol-v160*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_TPDL1_atezo_TPDL1) <- v160*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_CD28_CD80) <- v145*m_syn_T_APC*NAVG_nmol-v146*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_CD28_CD80_CD28) <- v146*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_CD28_CD86) <- v147*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_CD80_CTLA4) <- v148*m_syn_T_APC*NAVG_nmol-v149*m_syn_T_APC*NAVG_nmol-v151*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_CD80_CTLA4_CD80) <- v151*m_syn_T_APC*NAVG_nmol-v152*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_CTLA4_CD80_CTLA4) <- v149*m_syn_T_APC*NAVG_nmol-v150*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_CD80_CTLA4_CD80_CTLA4) <- v150*m_syn_T_APC*NAVG_nmol+v152*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_CD86_CTLA4) <- v153*m_syn_T_APC*NAVG_nmol-v154*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_CD86_CTLA4_CD86) <- v154*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_TPDL1_CD80) <- v157*m_syn_T_APC*NAVG_nmol-v158*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_TPDL1_CD80_TPDL1) <- v158*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_CD28) <- -v145*m_syn_T_APC*NAVG_nmol-v146*m_syn_T_APC*NAVG_nmol-v147*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_CTLA4) <- -v148*m_syn_T_APC*NAVG_nmol-v149*m_syn_T_APC*NAVG_nmol-v152*m_syn_T_APC*NAVG_nmol-v153*m_syn_T_APC*NAVG_nmol-v155*m_syn_T_APC*NAVG_nmol-v156*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_CD80) <- -v145*m_syn_T_APC*NAVG_nmol-v148*m_syn_T_APC*NAVG_nmol-v150*m_syn_T_APC*NAVG_nmol-v151*m_syn_T_APC*NAVG_nmol-v157*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_CD86) <- -v147*m_syn_T_APC*NAVG_nmol-v153*m_syn_T_APC*NAVG_nmol-v154*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_CTLA4_ipi) <- v155*m_syn_T_APC*NAVG_nmol-v156*m_syn_T_APC*NAVG_nmol
    d/dt(q_syn_T_APC_CTLA4_ipi_CTLA4) <- v156*m_syn_T_APC*NAVG_nmol
    d/dt(q_V_1_NabP) <- -v186-v187-v188
    d/dt(q_V_2_NabP) <- v186
    d/dt(q_V_3_NabP) <- v187

    # ---------- initial conditions (Supplementary Table S3 / rules) ----------
    q_V_C_nT0(0) <- 3706.89655172   # Table S3 (3706.9 cell)
    q_V_C_nT1(0) <- 2274.77477477   # Table S3 (2274.77 cell)
    q_V_P_nT0(0) <- 185344.827586   # Table S3 (185345 cell)
    q_V_P_nT1(0) <- 113738.738739   # Table S3 (113739 cell)
    q_V_P_Treg_CTLA4(0) <- (Treg_CTLA4_tot)*m_V_P*NAVG_nmol   # Table S6 rule 41
    q_V_T_C1(0) <- C1_0   # from initial_tumour_diameter via Eq 2 geometry; see Errata
    q_V_T_P0(0) <- 1e-09*vol_V_T0   # Table S3 (1e-18 molarity)
    q_V_T_P1(0) <- 1e-09*vol_V_T0   # Table S3 (1e-18 molarity)
    q_V_T_Treg_CTLA4(0) <- (Treg_CTLA4_tot)*vol_V_T0*NAVG_nmol   # Table S6 rule 40
    q_V_T_K(0) <- 27000   # Table S3 (27000 cell)
    q_V_LN_nT0(0) <- 155.689655172   # Table S3 (155.69 cell)
    q_V_LN_IL2(0) <- 1e-09*m_V_LN   # Table S3 (1e-18 molarity)
    q_V_LN_nT1(0) <- 113.738738739   # Table S3 (113.739 cell)
    q_V_e_P0(0) <- 1e-09*m_V_e   # Table S3 (1e-18 molarity)
    q_V_e_p0(0) <- 1e-09*m_V_e   # Table S3 (1e-18 molarity)
    q_V_e_P1(0) <- 1e-09*m_V_e   # Table S3 (1e-18 molarity)
    q_V_e_p1(0) <- 1e-09*m_V_e   # Table S3 (1e-18 molarity)
    q_A_e_M1(0) <- 0.0362959358945*m_A_e*NAVG_nmol   # Table S3 (2185.79 molecule/micrometer^2)
    q_A_e_M1p0(0) <- 1.66053906717e-11*m_A_e*NAVG_nmol   # Table S3 (1e-06 molecule/micrometer^2)
    q_A_e_M1p1(0) <- 1.66053906717e-11*m_A_e*NAVG_nmol   # Table S3 (1e-06 molecule/micrometer^2)
    q_A_s_M1(0) <- 1.66053906717e-11*m_A_s*NAVG_nmol   # Table S3 (1e-06 molecule/micrometer^2)
    q_A_s_M1p0(0) <- 1.66053906717e-11*m_A_s*NAVG_nmol   # Table S3 (1e-06 molecule/micrometer^2)
    q_A_s_M1p1(0) <- 1.66053906717e-11*m_A_s*NAVG_nmol   # Table S3 (1e-06 molecule/micrometer^2)
    q_syn_T_C1_PD1(0) <- ((T_PD1_total/A_Tcell))*m_syn_T_C1*NAVG_nmol   # Table S6 rule 16
    q_syn_T_C1_PDL1(0) <- ((C1_PDL1_base/A_cell))*m_syn_T_C1*NAVG_nmol   # Table S6 rule 20
    q_syn_T_C1_PDL2(0) <- (((C1_PDL1_base*r_PDL2C1)/A_cell))*m_syn_T_C1*NAVG_nmol   # Table S6 rule 21
    q_syn_T_C1_TPDL1(0) <- ((T_PDL1_total/A_Tcell))*m_syn_T_C1*NAVG_nmol   # Table S6 rule 19
    q_syn_T_C1_CD28(0) <- ((T_CD28_total/A_Tcell))*m_syn_T_C1*NAVG_nmol   # Table S6 rule 17
    q_syn_T_C1_CTLA4(0) <- ((T_CTLA4_syn/A_Tcell))*m_syn_T_C1*NAVG_nmol   # Table S6 rule 18
    q_syn_T_C1_CD80(0) <- ((C1_CD80_total/A_cell))*m_syn_T_C1*NAVG_nmol   # Table S6 rule 24
    q_syn_T_C1_CD86(0) <- ((C1_CD86_total/A_cell))*m_syn_T_C1*NAVG_nmol   # Table S6 rule 25
    q_syn_T_APC_PD1(0) <- ((T_PD1_total/A_Tcell))*m_syn_T_APC*NAVG_nmol   # Table S6 rule 28
    q_syn_T_APC_PDL1(0) <- ((APC_PDL1_base/A_APC))*m_syn_T_APC*NAVG_nmol   # Table S6 rule 32
    q_syn_T_APC_PDL2(0) <- (((APC_PDL1_base*r_PDL2APC)/A_APC))*m_syn_T_APC*NAVG_nmol   # Table S6 rule 33
    q_syn_T_APC_TPDL1(0) <- ((T_PDL1_total/A_Tcell))*m_syn_T_APC*NAVG_nmol   # Table S6 rule 31
    q_syn_T_APC_CD28(0) <- ((T_CD28_total/A_Tcell))*m_syn_T_APC*NAVG_nmol   # Table S6 rule 29
    q_syn_T_APC_CTLA4(0) <- ((T_CTLA4_syn/A_Tcell))*m_syn_T_APC*NAVG_nmol   # Table S6 rule 30
    q_syn_T_APC_CD80(0) <- ((APC_CD80_total/A_APC))*m_syn_T_APC*NAVG_nmol   # Table S6 rule 36
    q_syn_T_APC_CD86(0) <- ((APC_CD86_total/A_APC))*m_syn_T_APC*NAVG_nmol   # Table S6 rule 37

    # ---------- reported outputs ----------
    tumour_volume <- vol_V_T                                  # L
    tumour_volume_mL <- vol_V_T*1000
    tumour_diameter_cm <- 10*(6*vol_V_T/pi)^(1.0/3.0)         # dm -> cm
    cancer_cells <- C_total
    CD8_density <- x_V_T_T1/vol_V_T/1000                      # cell/mL
    CD4_density <- x_V_T_Th/vol_V_T/1000                      # cell/mL
    Treg_density <- x_V_T_T0/vol_V_T/1000                     # cell/mL
    MDSC_density <- x_V_T_MDSC/vol_V_T/1000                   # cell/mL
    PDL1_tumour <- x_syn_T_C1_PDL1_total
    Cc_atezo <- x_V_C_atezo                                   # nmol/L
    bsa_m2 <- BSA/100    # dm^2 -> m^2; nab-paclitaxel is dosed at 100 mg/m2
    Cc_nabp <- x_V_1_NabP                                     # ug/L
  })
}
