A multispecies framework for modeling adaptive immunity and immunotherapy in cancer

Predator-prey theory is commonly used to describe tumor growth in the presence of selective pressure from the adaptive immune system. These interactions are mediated by the tumor immunopeptidome (what the tumor “shows” the body) and the T-cell receptor (TCR) repertoire (how well the body “sees” cancer cells). The tumor immunopeptidome comprises neoantigens which can be gained and lost throughout tumorigenesis and treatment. Heterogeneity in the immunopeptidome is predictive of poor response to immunotherapy in some tumor types, suggesting that the TCR repertoire is unable to support a fully polyclonal response against every neoantigen. Importantly, while tumor and T-cell populations are known to compete with each other for intratumoral resources, whether between-lineage competition among peripheral T cells influences the TCR repertoire is unknown and difficult to interrogate experimentally. Computational models may offer a way to investigate these phenomena and deepen our understanding of the tumor-immune axis. Here, we construct a predator-prey-like model and calibrate it to preclinical and clinical data to describe tumor growth and immunopeptidome diversification. Simultaneously, we model the expansion of antigen-specific T-cell lineages and their consumption of both lineage-specific antigenic resources and lineage-agnostic, shared resources. This predator-prey-like framework accurately described clinically observed immunopeptidomes; recapitulated response-associated effects of immunotherapy, including immunoediting; and allowed exploration of treatment of tumors with varying growth and mutation rates.

Intratumoral predator-prey dynamics depend on the immunopeptidome. A. Simplified schematic of neoantigen presentation (created with BioRender.com). Neoantigens arriving from tumor-specific aberrations in DNA, inter alia, are presented on MHC molecules. CTLs with TCRs that recognize the neoantigen subsequently proliferate and attack the neoantigen-presenting cell. Not depicted: neoantigens arising from mechanisms other than somatic mutations (e.g., aberrant RNA splicing). pMHC, peptide-MHC complex. B. Schematic of immunodominance. Highly immunogenic neoantigens (dark green) that are also subclonal can be shielded from CTLs due to the abundance of more clonal neoantigens (light green). C. Schematic of interclonal competition. CTLs reactive against different antigens compete for the same pool of shared nonspecific stimulus, including cytokines and co-stimulation from antigen-presenting cells. DC, dendritic cell.
It is well known that some antigens are ignored by the immune system, often due to low abundance or presence in an immunologically privileged site. But it remains unclear why the TCR repertoire does not support a fully polyclonal response against every neoantigen during ICB treatment, clonal or not. We speculate this limitation could be due to the spatiotemporally restricted availability of resources and nutrients important for CTL survival and proliferation [16]. Cytokines critical to antigen-driven expansion (e.g., IL-2, IL-15), co-stimulatory signals from antigen-presenting cells (e.g., CD86 on dendritic cells), and other factors beyond direct antigen-mediated signaling through TCRs are a shared resource among CTL lineages (Fig 1C) [ [16][17][18]. The asymmetric and possibly contemporaneous presentation of multiple neoantigens might exhaust these nonspecific stimuli, leading to an inadequate immune response to some neoantigens. In humans, the clonal expansion of antigen-reactive CTL lineages can indeed suppress the expansion of other CTL lineages, a phenomenon called immunodominance that is well characterized within virology [19][20][21]. In cancer, ICB has been widely reported to affect the diversity of the CTL TCR repertoire, although whether it biases CTL expansion toward polyclonal or oligoclonal states appears to vary by ICB target and tumor histology [22]. However, given the shared nature of nonspecific stimuli, highly clonal tumors might be more susceptible to CTL-based immunotherapies due to reduced interclonal competition for nonspecific stimuli as compared to highly subclonal tumors rife with lower-quality, "distracting" neoantigens [16].
Neoantigen quantity is therefore critical to consider alongside immunogenicity [15,23]. However, the quantitative relationship between immunopeptidome diversity, CTL diversity, and response to ICB remains open to characterization. In this work, we construct a mathematical model of nonlinear predator-prey-like dynamics to describe the relationship between immunopeptidome diversity and CTL responses while accounting for interclonal CTL competition for nonspecific stimuli (Fig 2). Tumor growth and immunopeptidomics are modeled and calibrated with experimental data from two different species and three types of cancer, including an orthotopic murine model of bladder cancer, a murine model of colorectal cancer, a 100-patient cohort of non-small cell lung cancer (NSCLC) patients who underwent extensive immunopeptidome profiling [10,11,[24][25][26]. Simulated NSCLC-like immunopeptidomes and TCR repertoires are cross-validated on an independent cohort [27] and subject to a local sensitivity analysis. Finally, tumor response dynamics in the NSCLC-like cohort, as well as in two additional cohorts resembling other tumor types, are simulated under ICB treatment to evaluate pre-treatment biomarkers of response within the immunopeptidome and the TCR repertoire.

Model recapitulates tumor growth in immunocompetent and immunodeficient mice
BBN963 orthotopic bladder cancer cells generated in [24] were implanted in immunocompetent C57BL/6 and immunodeficient NSG mice [28] (Fig 3A). Thirty-four neoantigens putatively expressed by the cells were identified and tested for immunogenicity ex vivo [28] (see Methods). We used the results of these immunogenicity assays to scale in silico immunogenicity and assign one of two murine MHC alleles for presentation ( Fig 3B). Five-week simulations of tumor growth in the absence of CTL expansion and CTL-mediated killing accurately reproduced BBN963 growth dynamics in NSG mice (Fig 3C). Similarly, five-week simulations with CTL expansion and cytotoxicity restored accurately reproduced BBN963 growth dynamics in C57BL/6 mice ( Fig 3D).
To test model fidelity against an independent dataset, we used the growth dynamics of MC38 murine colon cancer cells in immunodeficient and immunocompetent mice reported in [26]. Assuming the same distribution of neoantigen immunogenicity as BBN963 cells and changing no parameters except growth rate, the model accurately captured MC38 tumor growth in both immunocompetent and immunodeficient mice (S2B- S2D Fig). These preliminary results suggested the generalizability of our predator-prey-like framework and motivated further simulation of human tumor-immune dynamics.

Model recapitulates the tumor immunopeptidome of NSCLC patients
We substituted human parameter values for steady-state T-cell abundance and lifespan from the mouse model, while also re-estimating tumor carrying capacity and T-cell sensitivity to antigenic stimulus. No other parameters related to tumor or T-cell dynamics were changed (Table 1). This enabled successful simulation of the growth and immunopeptidome of 100 human NSCLC tumors (Fig 4A) [25]. Tumor growth was simulated until target population sizes approximating the tumor diameters of TRACERx 100 NSCLC patients were reached. In contrast to the murine simulations, cancer cells in the human simulations were allowed to stochastically gain new neoantigens, lose expression of existing neoantigens, and lose MHC expression (see Methods, S2A Fig) [10,11]. Simulations recapitulated the distribution of clonal, subclonal, and total neoantigen burden within this cohort ( Fig 4B).

PLOS COMPUTATIONAL BIOLOGY
Interestingly, we observed substantial tumor and CTL diversification over time, as depicted in Fig 4C-D. The number of unique tumor lineages was consistently higher than the number of newly tumor-reactive CTL lineages due to lineage branching events caused by neoantigen and MHC loss. These cases represent the formation a new cancer lineage without introduction of a new neoantigen and its cognate predator CTL (e.g., expression loss of a current neoantigen). There was also a mildly positive correlation between the number of unique neoantigens and CTL diversity (Fig 4E). Murine tumor growth dynamics were recapitulated in divergent immune contexts. A. Schematic of NSG and C57BL/6 mouse inoculation with subcutaneous BBN963 tumors (created with BioRender.com). Tumor volumes were measured over five weeks [28]. B. Saito et al. identified putative neoantigens expressed by BBN963 and the MHC alleles predicted to present them (H2-Db, H2-Kb). Baseline-corrected ELISpot scores were used to scale neoantigen immunogenicity to the range [0,1]. C. Simulated tumor volumes (left, black) normalized to baseline volume (dashed gray line) compared against growth dynamics in NSG mice reported in [28] (gray circles and error bars). Nonspecific stimulus (center) and CTL populations (right) are shown. Dashed lines, 1. D. Simulated tumor volumes (left, black) normalized to baseline volume (dashed gray line) compared against growth dynamics in C57BL/6 mice reported in [28]. (gray circles and error bars). Nonspecific stimulus (center) and CTL populations (right) are shown. Black dashed lines, 1; Red dashed lines, approximate total abundance of 3 Ag-specific T-cell lineages observed in mice [28].

Simulated ICB highlights pretreatment clonal neoantigen burden as a biomarker of response
Prolonged antigen exposure results in CTL exhaustion, diminishing their proliferative and cytolytic capabilities. This process is accelerated by checkpoint molecules, many of which are overexpressed in the tumor microenvironment. ICB blocks these checkpoint molecules, reinvigorating the immune response and promoting the expansion of a select subset of CTL lineages ( Fig 5A) [29]. To simulate ICB of late-stage NSCLC, we constructed a 100-tumor virtual cohort comprising tumors of at least~5 cm in diameter. This roughly corresponds to the advanced stage of NSCLC for which ICB is indicated. Similar to the TRACERx 100 virtual cohort, we observed a mild positive correlation between neoantigen burden and CTL diversity prior to treatment ( Fig 5B). Various strengths of ICB treatment were explored within the model by simulating tumor growth for one year in the presence of 2-, 6-, and 10-fold increases in CTL sensitivity to antigen stimulus with concurrent 2-, 6-, and 10-fold increases in steady-state CTL carrying capacity. As expected, the proportion of patients with maximal tumor shrinkage of at least 30%-"responders"-increased with treatment strength [30]. Notably, response rates in the 10-fold treatment simulation approached levels observed for ICB in Phase 3 clinical trials ( Fig 5C) [31].
Pretreatment neoantigen burden correlated positively with deeper responses, in line with previous reports (Fig 5D) [5]. Interestingly, this was driven primarily by clonal neoantigen burden. No significant relationship was found between response depth and the number of unique subclonal neoantigens. We also examined whether the immunogenicity or clonality of the most immunogenic neoantigen (i.e., rank 1) could predict response. In this case, clonality outperformed immunogenicity in predictive sensitivity and specificity; however, both were outperformed by the product of the two, termed clonality-weighted immunogenicity ( Fig 5E). Responding tumors also had consistently lower Shannon entropy in their immunopeptidomes, corresponding to more homogeneous sets of neoantigens ( Fig 5F). These data support the importance of neoantigen clonality and highlight intratumoral heterogeneity as a key detractor from ICB response [5,32].

CTL dynamics correlate with ICB response and immunoediting
While responders tended to have greater pretreatment CTL diversity, further CTL diversification induced by ICB treatment strongly correlated with shallower responses (Fig 6A-B). Instead, expansion and maintenance of CTL lineages that were highly abundant before ICB initiation tended to elicit more profound responses ( Fig 6C) [22,33,34]. These findings support the notion that prognoses are most favorable when ICB induces focused, oligoclonal responses against highly clonal neoantigens, at least for certain tumor types [22,35].
Finally, we evaluated the potential for immunoediting against the rank 1 neoantigen during treatment. Immunoediting is the phenomenon by which cell lineages bearing certain neoantigens are selected against and killed by CTLs, resulting in diminishment or even eradication of certain neoantigens from the immunopeptidome [6]. Across the simulated tumors, we observed no correlation between the immunogenicity nor clonality-weighted immunogenicity of the rank 1 neoantigen and on-treatment immunoediting. There was, however, a very modest negative correlation between the immunoediting and the clonality of the rank 1 neoantigen alone (Fig 6D-6F). This suggests that clonality may supersede immunogenicity in provoking lineage-specific selection from CTLs; however, given the very modest effect and borderline significance, further work in larger patient populations-simulated or otherwiseis needed.

Model enables exploration of diverse tumor immunopeptidomes
To evaluate model generalizability, RCC-like and MSI-like immunopeptidomes with adjusted tumor growth and mutation rates were simulated and subjected to ICB (Fig 7A-7B) [32,36]. Pearson's ρ correlation between unique neoantigen number and tumor-reactive CTL Shannon entropy in the simulated 100-tumor cohort prior to treatment. Excludes one tumor with > 600 neoantigens for visualization purposes. C. Maximum tumor shrinkage under increasing levels of treatment strength: 2-(left), 6-(center), and 10-(center) fold increases in sensitivity to antigen stimulus and systemic carrying capacity were simulated. Each bar corresponds to a tumor. The dashed line at corresponds to the 30% or greater tumor shrinkage required to classify as partial response per RECIST 1.1 criteria [30]. D. Pearson's ρ correlation between response and pretreatment total, clonal, and subclonal neoantigen burden. Excludes one tumor with > 600 neoantigens for visualization purposes. E. Response status (colored bars) relative to rank 1 neoantigen immunogenicity (top), clonality (middle), and clonalityweighted immunogenicity (bottom). Sensitivity and specificity were calculated from the correctly classified fraction of responders (N = 36) and non-responders (N = 64), respectively. F. Fraction of responses among the X% most and least entropic immunopeptidomes, X being indicated by the X axis value. Entropy refers to Shannon entropy. Dashed lines, means.
https://doi.org/10.1371/journal.pcbi.1010976.g005 Spearman's ρ rank correlation between pre-treatment immunogenicity of the rank 1 most immunogenic neoantigen prior to treatment and the extent of immunoediting after one year of ICB treatment. Immunoediting was determined by calculating the fold change of the number of cells expressing the rank 1 neoantigen between the first and last day of ICB. E. Spearman's ρ rank correlation between pre-treatment clonality of the rank 1 most immunogenic neoantigen prior to treatment and the extent of immunoediting after one year of ICB treatment. F. Spearman's ρ rank correlation between pre-treatment clonality-weighted immunogenicity of the rank 1 most immunogenic neoantigen prior to treatment and the extent of immunoediting after one year of ICB treatment.

PLOS COMPUTATIONAL BIOLOGY
As expected, lower tumor growth and mutation rates led to lower neoantigen burden in RCClike tumors (Fig 7B). This resulted in poorer predicted responses in these tumors when compared to NSCLC-like and MSI-like tumors, in alignment with clinical experience (Fig 7C) [31,37,38]. Clonal neoantigen burden continued to correlate with response depth, although these results may have been skewed in RCC-like tumors by the high proportion of nonresponders ( Fig 7D).
Generally similar results were observed for neoantigen-rich MSI-like tumors. However, despite their higher neoantigen burden, the most heterogeneous MSI-like tumors were less likely to respond than most homogenous tumors ( Fig 7E). Furthermore, the CTL lineages that were most abundant prior to treatment in MSI-like tumors expanded less profoundly upon ICB treatment (Fig 7F). This is in sharp contrast to NSCLC-like tumors, in which oligoclonal CTL expansion was observed even in nonresponders. The nonlinear relationship between CTL expansion and neoantigen burden warrants further investigation.

Discussion
Characterizations of the quantitative relationship between the immunopeptidome and CTL clonal dynamics remain lacking. This work provides a mathematical, predator-prey-like framework capable of modeling the co-evolution of the tumor immunopeptidome and TCR repertoire. Built initially upon mouse data, this framework was straightforward to adapt for human tumor simulations. This allowed us to interrogate several clinically observed phenomena, including both pre-and on-treatment changes to the immunopeptidome and TCR repertoire. While several noteworthy computational models of neoantigen evolution have been published, these models consider either a fixed cancer cell population size or calculate death rates as a function of the cumulative immunogenicity of neoantigens expressed by a lineage or clone [2,15,39,40]. A recently published study modeled anti-tumor immunity against a growing tumor and incorporated both a dynamic immunopeptidome and neoantigen-specific CTL infiltrate [41]. No work, to our knowledge, has considered the possible interaction between tumor-reactive and tumor-ignorant CTL populations on a systemic level. Our results support the directive that neoantigen discovery methods should focus not only on immunogenicity, but also on clonality [42].
One key insight from this work is the shift away from neoantigen "dilution" as the predominant sequelae of intratumoral heterogeneity and toward CTL "distraction". Several studies report the ability of subclonally expressed neoantigens to dilute an otherwise strongly immunogenic neoantigen; that is, to shield it from immunoediting and dampen anti-tumor killing [8,14]. Recent modeling work suggests that highly mutable cancers can only grow and persist despite high neoantigen burden due to their ability to suppress anti-tumor immunity [43]. Furthermore, it was found that both greater subclonal neoantigen diversity and greater abundance of poorly immunogenic tumor lineages detract from anti-tumor immunity [41]. Based on these and our results, it is conceivable that a diverse immunopeptidome itself directly contributes to this suppression.
Notwithstanding the foregoing, our model supports the importance of neoantigen dilution and shows that neoantigen clonality may be more important than immunogenicity in predicting immunoediting and conceiving biomarkers of ICB responders. Considering clonality and immunogenicity in tandem is especially informative to ICB response (Fig 6E). However, while these metrics are useful, they largely ignore the CTL compartment. Holistically characterizing the anti-immunopeptidome response requires concomitantly assessing how the TCR repertoire evolves across space, time, and tumor types.

PLOS COMPUTATIONAL BIOLOGY
In situations where antigen stimulus is "shocking"-e.g., vaccination, or introduction of anti-CD19 CAR T cells into patients with B-cell cancers-clonal T-cell expansion can be as high as 3 logs over~2-4 weeks [44,45]. Conversely, chronic and gradual antigen stimulus like tumorigenesis does not necessarily provoke as harsh and rapid a T-cell response; otherwise, the tumor would be eliminated. Tumor microenvironments may also include substantial populations of exhausted T cells, which proliferate more slowly (or not at all) with antigen stimulus. To our knowledge, a dataset comprising patient TCR repertoires well before and soon after cancer diagnosis does not exist. This hinders our ability to quantify the effective T-cell expansion rate in humans upon gradual antigen stimulus in a clean and minimally confounded way.
However, pre-and post-anti-PD-1 therapy T-cell clonal expansion studies do exist [6]. We note that Anagnostou et al. 2017 [6] detected clonal expansion of the single most abundant T-cell clones in patients who responded to 10-15 months of anti-PD-1 therapy in the 25-to 44-fold range. Our results, which represent the top 5 most abundant T-cell clones in patients who responded to 12 months of anti-PD-1 therapy, closely approximate this, with peak expansion in the 15-to 20-fold range (Fig 6C). In addition, Anagnostou et al. 2017 [6] detected a decline in the frequency of these clones at the time of therapeutic resistance. We also found this to be the case in our simulations.
Notably, we did not differentiate between tumor-infiltrated and peripherally circulating T cells. Yost et al. previously demonstrated the concept of clonal replacement in solid tumors [46]. This is the idea that, under conditions of enhanced anti-tumor immunity like immunotherapy, extra-tumoral T cells replace (presumably already dysfunctional) intratumoral T cells to lyse tumor. This is supported by findings from Wu et al. 2020 [47]. We find this intriguing considering (1) the remarkable velocity of T cell trafficking naturally carried out for antigen surveillance, and (2) the similarly remarkable velocity of intratumoral T cell recirculation [48]. Put metaphorically, one might imagine the body as one contiguous lymphocyte highway, with tumors and peripheral tissues acting as "toll booths". At these toll booths, T cells stop but almost always proceed to re-enter the highway, except under rare and extenuating circumstances (e.g., recognition of cognate Ag).
Our model parameters (Table 1 and Eq 1) imply that tumor-reactive CTLs need to outnumber tumor-ignorant CTLs in order for tumor shrinkage to occur. Is this biologically plausible? This is challenging to unequivocally answer from published datasets, as most focus on TCR diversity indices instead of the quantification of bona fide, functionally validated, neoantigenreactive T cell frequencies. Some insight is provided by Wu et al. 2020 [47], who show that "dually expanded" (abundant in tumor and periphery) T cells can sometimes outnumber all other T cell clones combined. Of note, this was in treatment-naïve patients, for whom spontaneous tumor regression is presumably rare; one could imagine tumor-reactive clones further expanded under immunotherapy regimens. It seems plausible that the abundance threshold implied by our parameters is within biological ranges.
Here, we propose that competition for nonspecific stimulus by different CTL lineages is a mechanism for ICB resistance and provide a mathematical modeling framework for its investigation. Our simulations support interclonal CTL competition as a potential reason for why highly heterogeneous, subclonal tumors are associated with poor outcomes. We show that CTL diversification during treatment correlates with shallower responses,. In addition, many responding MSI-like tumors that responded to ICB had noticeably lower tumor-specific CTL expansion despite greater clonal neoantigen burden. Given these findings, "neoantigen dilution" may be a slight mischaracterization that diverts attention from the fact that CTLs, not cancer cells, effect anti-neoantigen immunity and tumor debulking.
Other limitations of this work include its spatially ignorant modeling of metastatic disease as one lumped tumor. Neoantigens can differ between patients and among metastases, potentially leading to increased systemic heterogeneity but decreased local heterogeneity [49]. Moreover, whether the immunopeptidome of individual metastases influences the already spatially biased trafficking of CTLs remains an open question. Are neoantigen-reactive CTLs more likely to linger in tumors where their target neoantigen is present? What if the neoantigen is highly subclonal-what is the probability they will physically encounter it? Is there a critical clonality threshold? These are critical questions that the model is not currently equipped to address.
Our work is primarily informed by data from two mouse models and two NSCLC studies. While the role of neoantigen clonality in NSCLC is generally consistent with its role in melanoma and colorectal cancer, it could differ in other tumor types such as low-grade glioblastoma and prostate cancer [32]. Moreover, the simulated response rate in MSI-like tumors was above that observed in the clinic, suggesting variability in treatment or CTL potency across tumor types. Ongoing clinical studies profiling the effect of ICB on TCR repertoire diversity in CRC may shed additional light on the veracity of our observations [50]. More generally, studies considering the coevolution of the cancer immunopeptidome and TCR repertoire diversity (or lack thereof) will be critical to achieving long-term disease control and cure.
In summary, this predator-prey-like framework accurately described clinically observed immunopeptidome diversity and recapitulated clinical evidence supporting the correlation between clonal neoantigen burden and ICB response. By including a dynamic CTL compartment, the model could also describe the effects of TCR diversity on predicted responses. Analysis of the tumor-immune axis also suggested that neoantigen immunoediting may correlate more strongly with pretreatment clonality than immunogenicity. This model can be readily adapted to future data published in other tumor types through simple, biologically motivated adjustments to parameter values. Further refinement could facilitate low-cost theoretical explorations of histology-specific therapeutic strategies in advance of developing and deploying them in the clinic. This work supports the further development of predator-prey-like models to elucidate nonlinear tumor-immune dynamics for clinical translation and utility.

Tumor growth dynamics
Tumors grew according to a modified logistic growth model, where N is the number of cells in the tumor, N i is the number of cells belonging to tumor lineage i, k g is the crude birth rate, K is the population carrying capacity, k d,i is the lineage-specific death rate (see Lineage-specific killing), and j is number of unique cancer cell lineages in the tumor: Modified logistic growth model describing tumor growth.  [51]. During mouse simulations, tumors were initiated from a single lineage~4.77 x 10 7 cells, corresponding to an approximate volume of 200 mm 3 at which tumor measurements began in [24,28]. During human simulations, tumors were initiated from a single founder cell with a gamma-distributed number of neoantigens [10]. Furthermore, the founder cell was assigned a 14% chance of having an MHC loss event [11] in one of its three MHC alleles. Further details on neoantigen immunogenicity, expression, and MHC presentation are provided in the Neoantigen features section. Parameter values, definitions, and their sources are provided in Table 1.

CTL clonal dynamics
At baseline, we assumed a tumor-ignorant and nonspecific CTL population equal to the steady-state population of T cells CTL ss in either mice or humans (Table 1). With the number of unique neoantigens k known from either a founder lineage (mice) or cell (humans), we then created k tumor-reactive CTL lineages (Fig 2). For simplicity, we assume TCRs are specific to a single neoantigen, which aligns with the high specificity of TCRs despite their considerable cross-reactive potential [52]. For each neoantigen among k neoantigens, a CTL lineage of 2 n cells was created, where n is a Poisson-distributed number n~P(μ). This represents the random number of divisions the founder CTL cell underwent prior to thymic efflux [53]. The remaining nonspecific CTLs were lumped into a tumor-ignorant, nonspecific population with a starting cell number equal to CTL ss minus the total number of neoantigenreactive CTLs.
Tumor-specific CTL lineages grew or shrank according to the quantity of their cognate antigenic stimulus and the availability of shared nonspecific stimuli, where CTL k is the kth CTL lineage reactive to the kth neoantigen, T d is its baseline turnover rate, C is the abundance of nonspecific stimulus, S k is the immunogenicity-weighted number of cancer cells expressing the kth neoantigen (see Neoantigen features), and T a is a scaling factor relating the rate of CTL k expansion per cancer cell expressing the kth neoantigen: Clonal dynamics of neoantigen-reactive CTL lineages.
This structure serves two purposes. The T d (C-1) term reflects clonal contraction when nonspecific stimulus is low and C < 1; depending on biological context, this could reflect the contraction of CTL clones upon antigenic threat elimination, as well as general immunodeficiency from low homeostatic cytokines. The S k T a C term reflects the influence of antigenic stimulus S k (scaled by T a and contemporaneously available co-stimuli C) on CTL expansion. The quantity of antigenic stimulus available to each CTL is assumed to decrease as the number of CTLs competing for a specific antigen increases hence the denominator of CTL k . Nonspecific CTLs presented with no antigen stimulus grew or shrank according only to the availability of shared nonspecific stimuli, where NSCTL is the number of non-specific CTLs: Growth dynamics of nonspecific CTLs.
Nonspecific stimulus C was consumed by T cells and non-T cells and replenished at rate k in .
The fraction of C consumed by cells other than CTLs was f out ; this was captured by the term k in f out . The remaining term k in (1 -f out ) representing the fraction of C consumed by CTLs was scaled by the number of CTLs at any given time relative to the number of CTLs at homeostasis, CTL ss . C therefore decreased when the total number of CTLs was greater than CTL ss and increased when it was less than CTL ss : Availability of nonspecific stimulus. ð4Þ Shannon entropy among the tumor-reactive CTL lineages was calculated using MATLAB's wentropy function.

Neoantigen features
Each neoantigen was assigned a non-negative real number R between 0 and 1 to describe its immunogenicity. In the mouse simulations, R for each of the 34 neoantigens identified in BBN963 cells in [24] was assigned by scaling the baseline-corrected ELISpot scores to the range [0, 1]. In the human simulations, R was assigned for each neoantigen in the founder cell by first sampling values for each of the TESLA criteria. Neoantigens were first assigned random values for MHC binding affinity, MHC binding stability, agretopicity, transcript abundance, and foreignness drawn from distributions defined by 146 NSCLC-derived peptides (S1 Fig) [4]. To achieve this, we fit an exponential distribution to transcript abundance and a zero-and oneinflated binomial distribution to foreignness. To account for correlation among MHC binding affinity, MHC binding stability, and agretopicity values, we fit a Gaussian copula to log-transformed values using MATLAB's copulafit function. Parameters associated with these fits are provided in Table 1. Transcript abundance was discarded due to conflation with clonality. The remaining four parameters were assessed against the following TESLA criteria: MHC binding affinity < 34 nM, MHC binding stability > 1.4 hours, agretopicity > 0.1, and foreignness < 10 −16 . Neoantigens that failed to satisfy these cutoffs were assigned R = 0, or no immunogenicity. Neoantigens that passed these cutoffs were assigned a value for R from an exponential distribution: Empirically determined neoantigen immunogenicity distribution.
1 � x � 0 In addition to intrinsic immunogenicity, each neoantigen was also assigned values to indicate its expression status and MHC binding partner. In all cases, neoantigens were assumed to be expressed at baseline. Murine simulations assumed each neoantigen could bind to one of two MHC alleles (H2-Db and H2-Kb), both of which were intact and expressed at baseline. Human simulations assumed neoantigen binding to one of three MHC alleles (HLA-A, HLA-B, and HLA-C), all of which were expressed at baseline, except in~14% of tumors where loss of one of three alleles at random conferred clonal MHC loss [11]. Thus, inclusion of a neoantigen in its parent cell's immunopeptidome depended upon both its own expression status as well as the intactness of its presenting MHC allele (S1 Fig). S k , the immunogenicity-weighted number of cancer cells with neoantigen k in their immunopeptidome, was used to stimulate the proliferation of CTL k per CTL clonal dynamics. It was calculated by multiplying the intrinsic immunogenicity of a neoantigen k by the number of cells whose immunopeptidomes included it, where E k is the lineage's expression status of neoantigen k, MHC k is the lineage's expression status of the MHC binding partner of neoantigen k, and R k is the intrinsic immunogenicity of neoantigen k: Neoantigen-driven CTL expansion.
In mouse simulations, cancer cell lineage immunopeptidomes were held constant within the short simulated timeframe due to the lack of data to inform event rates. In human simulations, stochastic modifications to a cancer cell lineage's immunopeptidome were allowed to occur. These events were one of three types: (1) neoantigen gain, (2) neoantigen loss, and (3) MHC allele loss. Neoantigen gain events, representing the generation of new and heritable neoantigens arising from DNA or other sources (e.g., aberrant RNA splicing) involved drawing TESLA parameters and assigning R immunogenicity values for each new neoantigen as described above. Neoantigen loss events, representing the genetic, epigenetic, transcriptomic, or metabolic suppression of a neoantigen, were represented as an irreversible change in the neoantigen's expression from 1 to 0; functionally, this removed the neoantigen from the lineage's immunopeptidome and reduced its S k to 0. MHC loss events, analogous to neoantigen loss events for MHC alleles, were represented as an irreversible change in one of the three MHC alleles' expression from 1 to 0. Functionally, this removed all neoantigens presented by that MHC allele from the lineage's immunopeptidome and reduced their S k to 0. Lineages that experienced one of these stochastic events were deducted 1 cell from their population size. An near-identical lineage with a modification corresponding to the stochastic event was then created with N = 1.
Neoantigen and MHC loss rates were directly calculated from the literature [10,11]. For example, only~33% of the neoantigens identified by Rosenthal et al. 2019 [10] were confirmed to be expressed, allowing estimation of a neoantigen loss rate from the relationship 0.67 = gain rate / loss rate. The MHC loss rate was set equal to the loss rate, assuming similar stochastic mechanisms would be responsible for the loss of either from the transcriptome or proteome. The gain rate was manually estimated by scaling until the distribution of subclonal neoantigen repertoires for our virtual TRACERx 100 cohort were comparable to the clinically observed data. Clonal neoantigen burden was directly assigned at t = 0 from the distribution of clonal neoantigen counts at baseline in the TRACERx 100 cohort; no fitting was performed. Finally, the clonal and subclonal neoantigen quantities were summed and compared to clinical data on a total neoantigen basis.

Lineage-specific killing
Each unique clone in the tumor was killed at an independent and saturable rate calculated from the contents of its immunopeptidome and the abundance of cognate, reactive CTLs.
Here, k d,i is a dimensionless value used to modify the logistic growth described in Eq 1; f die is the fraction of cell death not attributed to CTL killing, k mag and k kill are two parameters defining the saturable killing function. k mag determines whether k d,i can exceed 1, resulting in tumor shrinkage. j is the number of unique neoantigens harbored by the lineage: Empirically derived function for cancer cell lineage-specific death rates.
; * k kill ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where ; ¼ There two components of these death rates are natural turnover f die on the left and antigenspecific killing on the right. Antigen-specific killing is driven by ;, which itself is calculated from the quantity of CTLs reactive against the j neoantigens on a tumor lineage (as a fraction of total CTLs). CTL k is only counted if the kth neoantigen and its MHC are expressed (E k and MHC k = 1). k mag controls the maximum velocity of antigen-specific killing. Assuming k mag = 1, i.e., complete eradication of the tumor within one day, k kill modulates how rapidly k d,i approaches k mag with increasing values of ;.

Simulated tumor growth in mice
BBN963 orthotopic bladder cancer cells were previously established by chronic exposure of C57BL/6 mice to 0.05% N-Butyl-N-(4-hydroxybutyl) nitrosamine [24]. Neoantigens were predicted computationally and screened for bona fide immunogenicity in an ex vivo ELISpot assay (S2A Fig). Additional details can be found in [24]. BBN963 cells were implanted in immunocompetent C57BL/6 or immunodeficient NSG mice [28]. Tumor volumes were measured weekly for five weeks once the tumors reached 200 mm 3 . Tumor growth volumes were digitized using WebPlotDigitizer [54] and used to calibrate mouse-specific model parameters in Table 1. An independent dataset describing the growth dynamics of MC38 murine colon cancer cells in immunodeficient and immunocompetent mice was used to validate the mouse model. MC38 tumors were subcutaneously established in immunocompetent C57BL/6 or immunodeficient RAG1 KO mice [26]. Whole exome and RNA sequencing were used to identify 489 expressed neoantigens. In the absence of reported MHC allele specificity and ex vivo-validated immunogenicity, each neoantigen was assigned a random for MHC allele specificity and immunogenicity from the distribution reported for BBN963 cells [24]. Tumor growth in both immunodeficient and immunocompetent mice was successfully recapitulated by using the same parameters as BBN963 cells, except growth rate, which was increased 50% (S2B-S2D Fig).

Simulated tumor growth in humans
The TRACERx 100 cohort [25] contained 100 patients with varying stages of NSCLC: 26 stage I, 36 stage IB, 13 stage IIA, and 25 stage IIB/IIIA/IIIB. These patients underwent multiregion biopsy for extensive immunopeptidome profiling [10,11]. Tumor diameters were estimated for each stage using conventional NSCLC staging criteria [55]. Clonal and subclonal neoantigen burden digitized from [10] with WebPlotDigitizer, as well as reported rates of neoantigen expression loss and MHC clonal/subclonal loss of heterozygosity, were used to calibrate simulated immunopeptidomes. In total, 100 tumors were simulated to target diameters mirroring the TRACERx 100 cohort. Clonal and subclonal neoantigen burden for each tumor was then calculated by determining the number of neoantigens present in 100% or < 100% of lineagespecific immunopeptidomes.

Human immunopeptidome cross-validation and local sensitivity analysis
Simulated immunopeptidomes were evaluated against an independent cohort of 12 patients [27] (S3 Fig). Putative neoantigens from single-region biopsies collected from 12 patients prior to ICB treatment were assessed. Total neoantigens were calculated by filtering for positive netCTL classification. To apply the TESLA criteria, we calculated agretopicity by taking the ratio of somatic and wild-type MHC binding affinity. For predictions containing somatic peptides of identical sequences but different predicted MHC binding alleles, the more strongly binding prediction was used. Binding stabilities were predicted with NetMHCStab-Pan using default parameters [56]. The number of neoantigens that passed three of five the TESLA criteria (MHC binding affinity < 34 nM, MHC binding stability > 1.4 hours, and agretopicity > 0.1) was multiplied, in the absence of reported foreignness, by a correction factor equal to the ratio of neoantigens in the TRACERx 100 [10] cohort passing the MHC binding affinity, stability, and agretopicity criteria that also did or did not pass the foreignness cutoff of < 10 −16 . This reduced the number of predicted neoantigens in the immunopeptidome by 80.6%. Neoantigen clonality was not assessed due to potential sampling biases intrinsic to single-region biopsies. For each of the seven discrete tumor stages, tumor diameters were estimated for each stage using conventional NSCLC staging criteria [55]. A cohort of 40 tumors was simulated for each of the seven diameters. The simulated and experimentally reported immunopeptidomes are shown in S3

Exploratory simulations of ICB treatment and other tumor histologies
A cohort of 100 tumors was simulated until a target population roughly equivalent to a~5 cm diameter tumor. Treatment was implemented to increase the sensitivity of CTLs to antigenic stimuli as well as the population-level carrying capacity CTL ss . We assumed these two effects were of both driven by ICB.
Implementation of ICB within the model.
To emulate renal cell carcinoma (RCC)-like tumors, a cohort of 100 tumors with 50% lower growth rate and 90% lower neoantigen gain, neoantigen loss, and MHC loss rates was simulated. To emulate microsatellite instability-high colorectal cancer (MSI)-like tumors, a cohort of 100 tumors with 100% higher growth rate and 200% higher clonal neoantigen burden was simulated. ICB treatment was simulated in these tumors as they were for NSCLC, described in Eq 8.

Statistical analyses
All simulations and statistical analyses were performed in MATLAB R2020b. Pearson's and Spearman's ρ were used as appropriate to assess linear correlation among continuous or ranked variables. Continuous variables in the local sensitivity were assessed for statistical significance at α = 0.05 using one-way ANOVA and a post-hoc pairwise Student's t-test, with Bonferroni correction for multiplicity. Details of specific statistical analyses are provided in their corresponding figure captions.
Supporting information S1 Fig. Model-simulated neoantigens account for stochastic events and TESLA criteria. A. Each neoantigen was assigned values for the TESLA parameters drawn from the distributions for NSCLC tumors in [4]. TPM, transcripts per million. Distributions (top) and values (bottom) of sampled and experimental parameters from 146 NSCLC-derived peptides are shown. B. Schematic of stochastic neoantigen events (created with BioRender.com). A founder cell was established with a random number of neoantigens and assigned a neoantigen gain rate from a lognormal distribution. Each neoantigen was also assigned to be presented by one of three MHC alleles. Neoantigen and MHC expression could also be lost by a cell during tumor growth. (TIF)

S2 Fig. Characterization of the BBN963 immunopeptidome and validation in MC38 cells.
A. Schematic of BBN963 cells establishment by chronic exposure of C57BL/6 mice to 0.05% N-Butyl-N-(4-hydroxybutyl) nitrosamine (created with BioRender.com). Neoantigens were called and screened for bona fide immunogenicity in an ex vivo ELISpot assay. Additional details can be found in [24]. B. Growth dynamics of MC38 were simulated using the same parameters as BBN963 cells, except growth rate, which was increased 50% [26]. C. Simulated tumor volumes (left, black) normalized to baseline volume (dashed gray line) compared against growth dynamics in RAG1 KO mice reported in [26] (gray circles and error bars). Nonspecific stimulus (center) and CTL populations (right) are shown. Dashed lines, 1. D. Simulated tumor volumes (left, black) normalized to baseline volume (dashed gray line) compared against growth dynamics in C57BL/6 mice reported in [26] (gray circles and error bars). Nonspecific stimulus (center) and CTL populations (right) are shown. Dashed lines, 1. (TIF)

S3 Fig. Validation of simulated neoantigen burdens in an independent cohort.
A. Neoantigen predictions reported in [27] of 12 patients with sufficient pretreatment tissue to characterize candidate neoantigens were assessed. Total neoantigens were calculated by filtering positive netCTL classification. To apply Wells' criteria [4], we calculated agretopicity by taking the ratio of somatic and wild-type MHC binding affinity. For predictions containing somatic peptides of identical sequences but different HLA alleles, the stronger predicted binding was used. Binding stabilities were calculated with NetMHCStabPan using default parameters. As biopsies were from a single region, neoantigen clonality was not assessed. Tumor diameters were estimated for each tumor using conventional NSCLC staging criteria. For each of seven discrete estimated tumor diameters, a cohort of 40 tumors was simulated. (TIF) S4 Fig. Local sensitivity analysis of model parameters. A. Mean number of unique neoantigens relative to baseline (N = 100 simulations per condition). No significant fold differences in neoantigen number were found from pairwise Student's t-tests with nominal and Bonferroni-corrected significance thresholds. B. Mean TCR entropy relative to baseline (N = 100 simulations per condition). No significant fold differences in Shannon entropy were found from pairwise Student's t-tests with nominal and Bonferroni-corrected significance thresholds. (TIF)