Cancer-immune coevolution dictated by antigenic mutation 1 accumulation

Abstract


Introduction
The immune system does not only protect our body from infectious diseases caused by various pathogens, but also gives the first response against threats emerging within the body, such as cancer.Immune cells try to identify and eliminate tumour cells, which may express antigens not found on normal cells.Meanwhile, tumour cells attempt to hide or evade from immune surveillance [1].From an ecological perspective, this antagonistic relationship gives rise to complex dynamics between a tumour and its microenvironment [2,3].On an evolutionary level, these two cells types, though belonging to the same organism, coevolve, each adapting to genotypic and phenotypic changes in the other population [4].We are interested in the eco-evolutionary dynamics of genotype-specific interactions between cancer and immune cells arising from a continuous introduction of new antigens with the context of immune surveillance and escape.Dunn et al. described the battle between the immune system and emergent tumours in three stages, termed the three Es of cancer immunoediting [5].The first, Elimination, was formulated by Burnet as the immunosurveillance hypothesis, which states that the immune system can win this battle and eliminate small cancers [6,7].
According to the "bad luck" hypothesis [8], this is a frequent occurrence: only stochastically does the immune system allow cancers to sneak through [4].Should the cancer do so, it enters Equilibrium.Indeed, some cancers take years to grow to detectable size [5], and there is evidence for small persistent tumours coevolving with the immune system [9].Finally, Escape, when the cancer evolves mechanisms of evasion, and grows to a size that can be detected [10].Deciphering genetic footprints from this coevolutionary process is of utmost clinical relevance from prognosis to treatment [11], especially since current measures such as immune infiltration, evidence of immune escape and tumour mutational burden (TMB) are not foolproof markers of immunotherapeutic success and overall survival [12].
Mutations accumulated during cancer evolution increase intratumour heterogeneity, providing a wide landscape of genotypes to improve the tumour's persistence [13].Propensities for mutating more rapidly, increasing growth rate, developing metastases, and acquiring resistance to eventual treatment are possible consequences of this accrual of mutations [1].However, mutations may also alert the immune system of the presence of tumour cells and initiate a suppressive response.This occurs when a mutation acquired by a cancer cell (called an antigenic mutation) results in the presentation by human leukocyte antigen (HLA) of immunogenic peptides called neoantigens at the cell surface [14].These neoantigens are recognisable by cytotoxic T lymphocytes (CTLs), specialised effector cells of the adaptive immune system, which bind to the presented neoantigen via a T cell receptor at their surface and kill the targeted cancer cell [15].While antigenic mutations are neutral in the absence of an immune response [16], their fitness in general depends on the likelihood of neoantigen presentation and recognition by effector cells [17].
Antigenic mutations so targeted by the immune system are thus under negative selective pressure [18].The killing of cancer cells carrying antigenic mutations releases further neoantigens, leading to a cyclic process of generating immunity to cancer [19].Cancer cells, however, may in turn combat the immune response via several mechanisms.Cancer cells can inhibit effector cells, such as by expression of programmed cell death-ligand 1, which normally presents in healthy cells to stop the attack of immune cells and can exhaust CTLs interacting with cancer cells that carry it [20,21].Some cancer cells may escape the immune response, such as by reducing neoantigen presentation through down-regulating HLA [14] or by immune-editing and losing antigenic mutations due to the immune negative selection [22].Cancer cells can even develop immune exclusion, physically restricting the immune cells' access to the tumour [1].These processes may arise individually or in concert, and they mean the immune system plays a crucial role in shaping a tumour's evolution [1,10].
Immunotherapies leverage the immune system to reverse the evolution of evasion by mitigating or interfering with these processes [21].Cytokine-based and tumour-infiltrating leukocyte-based immunotherapies help increase the effector cell population [15,23]; immune checkpoint blockade therapies restore the effector response to immuneescaped cancer cells [12,19]; other immunotherapies simply boost the ability of CTLs to kill tumour cells [15,24].Understanding of the coevolution between cancer and immune cells, along with the cell-to-cell interactions that drive it, may inform when therapies will succeed-and why.
There has been a long history of studying antagonistic coevolution experimentally and theoretically [25][26][27][28], with extensive literature on mathematical models of tumour-immune interactions [21,29,30].Deterministic models often describe the antagonistic relationship between cancer and immune cells as obeying Lotka-Volterra dynamics, with the immune system predating on its tumour prey [31,32].Stochastic modelling of the tumourimmune system has been explored by e.g.George and Levine, who characterised escape as a random process with sequential mutations [4,33] and subsequently framed cancer evolution as an active optimisation process in response to an evolving immune landscape [34,35].In concert with patient sequencing data, Lakatos and colleagues proposed and applied a model of random antigenic mutation accumulation in order to describe the negative selection that neoantigens undergo, resulting in neutral-like evolutionary dynamics [10,11].Recently, Chen & Xie et al. included negative frequency dependence in a similar model of antigenic mutation accumulation, so that tumours are only immunogenic when a sufficiently-large proportion of their cells present neoantigens, predicting that tumours undergoing this frequency-dependent selection have poorer treatment outcomes than their negative selection counterparts [36].
Coevolution between the immune system and the threats it faces has often been studied in a gene-for-gene framework, which centres on the genetic makeup of multiple populations being tracked [37].However, there is a dearth of stochastic models that explore the explicit evolutionary dynamics of both tumour and immune cell populations [34]: the aforementioned models either encapsulate the immune response into a selection parameter s, implicitly assuming that effector cells react perfectly and instantaneously to a cancerous threat [11,36]; or, they are deterministic and thus miss out the critical impact of random processes on small population sizes, while omitting direct genetic information with which sequencing data can be compared [38].
We address this gap in the literature by presenting and analysing a novel stochastic coevolutionary model of tumour-immune dynamics.As in Lakatos et al.'s model, antigenic mutations accrue in cancer cells undergoing a branching process and are negatively selected against by the immune system [11].The adaptive immune system is represented by specialised effector populations that react to emergent neoantigens [39], leading to complex dynamics on both tumour and immune fronts.We focus on interactions between effector cells and their targets, incorporating both active and passive recruitment of effector cells, killing of cancer cells, and inhibition of effectors by cancer cells, while including explicit mechanisms of escape.In particular, we inspect how these interactions are central in determining the evolution of the system.Stochastic simulations allow us to characterise the various dynamics and outcomes that emerge, along with informing the impact of immunotherapy.

Birth-death process with mutation accumulation and escape
Consider a population of cancer cells undergoing a stochastic branching process with per-capita birth rate b and per-capita death rate d (see Figure 1A).We assume births and deaths happen randomly based on their rates, with exponentially-distributed wait times between events, as in the Gillespie algorithm [40].At each cancer cell division, daughter cells inherit their mother cell's mutations and acquire a random number of new mutations, drawn from a Poisson distribution with mean λ [41], as shown in Figure 1B.Similar to Lakatos et al., we focus on the exomic region [11], where the value of λ ranges between 1 and 10.All mutations are unique, as in the infinite sites approximation, where the exome is considered long enough for two co-occurring point mutations to be negligible [42].We distinguish between two types of mutations: with probability p a , a mutation is antigenic and can be recognised by given effector cells, and with probability 1 − p a it is neutral [11,15].During a division, there is a probability p e that an antigenic mutation escapes the immune system, therefore making the cell possessing it (and all of its descendants [36,43]) undetectable to the immune system (see Figure 1B).d, respectively (here ∅ represents no cell).B During a division event, daughter cells inherit all of their mother cell's antigenic and neutral mutations, depicted by numbers (where underlined numbers are antigenic).Cells carrying antigenic mutations have probability p e to escape and become neutral, as shown in the lower daughter cell with mutation 2. Each daughter cell also acquires a random number (drawn from a Poisson distribution with mean λ) of new mutations, where each mutation is antigenic with probability p a and neutral with probability 1 − p a .C For each antigenic mutation i present in the system, a corresponding effector cell population E i exists (blue), which grows with constant rate B and shrinks with per-capita death rate D. D Antigenic mutations in cancer cells (such as i and k) display unique neoantigens at the cell surface, whereas neutral mutations (such as j) do not.The neoantigens can be identified by specialised effector cells, which can only interact with the corresponding cancer cells.E When a cancer cell carrying the antigenic mutation i meets an effector cell of type i, three outcomes are possible: active recruitment of another effector cell of type i with rate α i , killing of the cancer cell with rate β i and inhibition/exhaustion of the effector cell with rate γ i .The expressions for the rates are found in equation (1).

Specialised effector response and cancer-immune interactions
Neoantigens trigger responses from the adaptive immune system: each antigenic mutation i calls a unique, specialised effector population E i , as in gene-for-gene coevolution [37].Whenever mutation i exists within the population, effector cells of type i are passively recruited from the body at constant rate B (∅ → (E i ), for ∅ the absence of a cell and (E i ) an effector cell of type i, with rate B), and die with per-capita rate D [44,45], as shown in Figure 1C.
To each antigenic mutation i, we associate two random numbers, each drawn independently from an exponential distribution with mean 1: an antigenicity A i , describing the propensity for an effector cell (E i ) to kill a cancer cell possessing mutation i, and an immunogenicity I i , which relates to the rate at which active recruitment of more effector cells of type i occurs during the specified tumour-immune interaction (see equation (1) and Figure 1E).
These can be thought to encapsulate the probability of an antigenic mutation leading to the presentation of neoantigens by HLA molecules and subsequent recognition by CTLs [17].
Correspondingly, when a cancer cell with an antigenic mutation presents a neoantigen at its surface, only effector cells of the corresponding population can interact with it [21], as shown in Figure 1D.Interactions between cancer cells possessing an antigenic mutation i and corresponding effector cells of type i have three possible outcomes: active recruitment with rate α i , killing with rate β i and inhibition/exhaustion with rate γ i (see Figure 1E): with rates where α 0 , β 0 , γ 0 and h α are constants for the entire tumour.Note that the reactions of active recruitment and inhibition/exhaustion have opposite effects on the effector cell population; here, we assume α 0 > γ 0 , the opposite of which is rarely considered [46], as it leads to net decreasing effector populations upon interactions with cancer cells.However, active recruitment α i obeys a type-II functional response rather than a type-I (that is, linear) response as inhibition/exhaustion γ i [47,48].This implies that (α i − γ i )C i is not linear in the cancer cell population C i , whereas if α i too satisfied a type-I response it would be.Instead, as shown in Figure A1 of the Appendix, (α i − γ i )C i is positive for small populations C i and negative once the clone C i is large enough.

Mutational distributions
We are interested in genetic information relating to each of the two populations, which we will use to identify and measure the coevolution between effector and cancer cells.For the cancer cells, we define the following summary statistics: the site frequency S j , denoting the number of antigenic mutations present in j cells, and the single-cell mutational burden B k , the number of cells that possess k antigenic mutations.The conglomeration of these (for non-negative integers j and k, respectively) form the site frequency spectrum (SFS) and the singlecell mutational burden distribution (MBD).These satisfy , where C denotes the total population of cancer cells and M denotes the total number of antigenic mutations [49].We will write U for this quantity (that is, for either side of the previous equality), the total number of antigenic mutational occurrences.
For the effector population, we define the average antigenicity ⟨A⟩ and the average immunogenicity ⟨I⟩ as follows: Models that focus only on the cancer cell populations ascribe an antigenicity to the tumour itself [11,36]; here, by considering both cancer and effector populations, we can also average the antigenicities and immunogenicities over all effector cells.By choosing to average over the effector cells as in equation ( 2), ⟨A⟩ and ⟨I⟩ become measures of the gene-for-gene immune response and its potency in fighting its tumour target, providing an additional angle on the coevolutionary system not present in single-population models.
We implement our model using a standard Gillespie algorithm [40].A realisation ends when either the cancer cell population exceeds a threshold size of K (which we call no suppression and write T K for the time that this occurs); or, the cancer cell remains below K (which we call suppression and write T end as the end of the Population of effector cells of type i, with total effector cell population E B, D 0.2, 0.1 Passive recruitment and per-capita death of an effector cell α i , β i , γ i -, -, 10 −3 Active recruitment, killing and inhibition/exhaustion rates, respectively, for interactions between cancer and effector cells of type i A i , I i Antigenicity and immunogeniticy, respectively, of antigenic mutation i (each drawn from Exp(1)), with averages over the effector population ⟨A⟩ and ⟨I⟩ S j Site frequency: number of antigenic mutations occurring j times B k Single-cell mutational burden: number of cancer cells with k antigenic mutations M Total number of antigenic mutations U Total antigenic mutational occurrences: Ending cancer population size for stochastic simulations T end 16 Maximum time for a simulation; if the cancer cell population never grows to size K, then the outcome is suppression, otherwise it is no suppression (at time T K ) Table 1.Notation and baseline parameter values, which are used for all simulations unless otherwise specified.

Cyclic dynamics between antigenicity and immunogenicity
Antagonism between effector cells and their cancer targets result in a range of complex dynamics.For elevated immune efficacy such as high values of α 0 and β 0 (the active recruitment and killing parameters), emerging tumours are suppressed in early stages.However, when the probability of escape p e is large enough, then neutrallike dynamics ensue, as most cancer cells evade the immune response.While it is expected that negative selection of the immune system on new antigens will selectively prune cancer cells with more antigenic mutations [11], this process depends on the effector populations themselves, through the interactions formulated in equation (1).
Since here we explicitly model both tumour and immune cell types, we are able to explore their population dynamics in concert.
Indeed, we can clearly observe a dynamical response between the cancer and effector cells in single evolutionary trajectories.Figure 2 depicts two representative realisations of our stochastic simulations, one for a lower mutation rate (λ = 1) and one for a high mutation rate (λ = 10), the latter of which can be thought of as representing hyper-mutated tumours [11].The population dynamics of Figures 2A and 2C show the effector population undergoing irregular spikes, out of phase with the cancer population, as seen in other coevolutionary models [50,51].Because of the specialised nature of our model, these arise when a single antigenic mutation i causes the rapid active recruitment of the corresponding effector cell population E i to combat the subclone possessing mutation i during that time period.Once mutation i is eradicated from the cancer cell population, the effector cells E i specialised to their neoantigen are removed from the system, as they no longer play a dynamical role and die out exponentially (see Figure 1C).The expected frequency and amplitude of these effector spikes can be approximated, as described in Section A8 of the Appendix.
We validate this by inspecting the coevolution of antigenicity and immunogenicity in the corresponding single realisations.In Figures 2B and 2D, we can see that the effector spikes arise due to the emergence of a mutation i that has a large immunogenicity I i (significantly greater than the mean ⟨I⟩), as demonstrated in equation (A18) of the Appendix.Unexpectedly, in addition to this elevated immunogenicity, Figures 2B and 2D imply that these spikes also arise for mutations with low antigenicity A i < ⟨A⟩.These fluctuations in ⟨I⟩ and ⟨A⟩ are coupled to the spikes in the population size; had a mutation emerged with a large antigenicity, it would have been quickly eradicated, and so the effector population size would not have grown sufficiently to be identified as a spike.As expected, the higher mutation rate (λ = 10) leads to more cycles of spikes in the cancer-effector dynamics as well as in the coevolution of antigenicity and immunogenicity dynamics.

Interactions dictate outcome of tumour progress
The outcome of the coevolution between the immune system and a tumour depends strongly on the interaction parameters between effector and cancer cells (equation ( 1)).Figures 3A and 3D illustrate heat maps for the tumour suppression proportion across different values of the active effector recruitment rate α 0 and the cancer killing rate β 0 .Due to the stochastic nature of the model, the outcome for a given parameter set is probabilistic.
As expected, for higher values of α 0 and β 0 , more effector cells are recruited and more cancer cells are killed, thus the suppression increased.While this pattern holds for different mutation rates (Figure 3), it is more dominating when mutation rate is high.For a higher mutation rate (λ = 10), the immune system is more effective since there are more antigenic mutations to target, resulting in a high suppression (see Figure 3D).
For a lower mutation rate (λ = 1), however, only a small proportion of tumours are suppressed even with strong immune response (Figure 3A).This is because only a fraction of the population is antigenic (and thus undergoing immunoselection): observe Figures A6A-C versus A6D-F of the Appendix.It also worth noting that in this case most suppressed tumours go extinct at early times (Figure A5A) due to stochasticity.
When the active effector recruitment rate is smaller than the effector inhibiting rate α 0 < γ 0 (Figure 3 red dashed line), the net outcome for effectors during an interaction is negative (though there are subtleties here due to the difference in functional responses; see Section A4 of the Appendix for details).This implies that the resulting suppressed tumours in this parameter region were either defeated by a passively-recruited effector population, or by effector types i corresponding to antigenic mutations with particularly high immunogenicities I i .This is possible since I i is drawn from an exponential distribution (and thus can be high), meaning that even when α 0 is low, the active recruitment can be greater than the inhibition/exhaustion by chance.This is exactly what we observed in our stochastic simulations.Under high mutation rate (λ = 10), most tumours are not suppressed under low interaction parameters α 0 and β 0 (Figure 3D), whereas for high interaction parameters most tumours go extinct (Figure 3E) rather than maintaining a slow growing pace (Figure 3F).Again, we see a similar pattern under a low mutation rate (λ = 1), though more weakly and with more noise.
Interestingly, for a high mutation rate there exists an intermediate range of interaction parameters that allows for tumours to neither go extinct nor to grow to capacity (Figure 3F).One further interpretation of Figure 3D is the presence of a killing threshold (here, at β 0 ≈ 10 −1 ) above which tumours are suppressed, no matter the active recruitment rate α 0 .The rest of the domain of Figure 3D then exhibits a much stronger dependence on the active recruitment rate α 0 , in line with the results of Wilkie & Hahnfeldt [16].), the rate of inhibition/exhaustion (see Section A4 for further discussion); see the Discussion section for details of point e.All parameter values not specified here are listed in Table 1.

Genetic markers of selection
Stochastic mutation accumulation in an exponentially growing population has been widely studied [52,53].
When the immune system has little impact on the cancer cell population, therefore, it is unsurprising to see consistent increases in the average number of mutations per cancer cell, as in Figure 4A.The theoretical expectation of neutrally accumulated mutations can be approximated by the average flux of number of new mutations entering the system, 2bt [54], as plotted.(Note that a correction for large times has recently been shown to apply [55], which helps explain why the dashed line over-estimates the simulated data in Figure 4.) When the effector population is selectively killing, however, possessing more antigenic mutations makes a cancer cell more likely to get killed, so the average number of antigenic mutations per cancer cell does not continuously increase with the population, as shown in Figure 4B.This effect of selection is even more pronounced in the case of a higher mutation rate (λ = 10), as shown in Figure A7   The average number of antigenic mutations per cell is the mean of the MBD: ⟨B⟩ = U/C, which, like the SFS, can be extracted from bulk data.However, single-cell sequencing data can also provide information on the MBD overall, which can be used in combination with bulk information to infer what selection is taking place in the system [56].As discussed in Section A5 of the Appendix, the expected neutral distributions of the SFS and the MBD have been solved [49,57], so divergence from these may demonstrate the presence of selection and the strength of cancer-immune interactions.In Figure 5, the SFS and MBD averaged over 100 realisations are plotted in conjunction with the corresponding theoretical predictions (black dashed line) for the case with no immune response (that is, where all mutations are neutral).The first and third rows (green points) represent neutral mutations, whereas the second and fourth rows (orange points) measure antigenic mutations; for the MBDs, the mean (that is, the value U/C) was plotted in dashed vertical lines for both the simulated data (in green and orange) and the theoretical predictions (in grey), which were calculated with equations (A12) and (A13) of the Appendix.
We notice that in the case of low immune effectiveness (Figures 5A-D), there is little deviation from the neutral expectation.When the immune system plays a larger role, however, the distinction is significant, as in Figures 5E-H.In particular, the cells with more antigenic mutations were more selectively killed by the immune system, so that the tail of the MBD is smaller and the mean is shifted down, as seen most prominently in Figure 5H.As before, when the mutation rate is higher (λ = 10), these effects are more striking, as shown in Figure A8 of the Appendix.On the other hand, the SFS shows limited difference from its neutral expectation (Figures 5A and 5E), reiterating the importance of integrating single-cell data with bulk sequencing data in identifying immune effects.Under stronger negative selection, we expect a depletion of antigenic mutations from the population and thus lower site frequencies than the theoretical prediction [11].This is visible as a slight depression of the data compared to the neutral prediction in Figure 5F, though such an observation is much clearer for a faster immune response, as discussed in Section A9.  1.

Discussion
Coevolution between effector cells and their cancer targets sets the stage for the emergence and subsequent development of a tumour.Based on expansive literature in this field [11,31,33], we focus on important perspectives which have not yet been addressed by the previous models: in particular, the stochastic nature of these early-stage small cell populations as well as explicit interactions between cancer cells and immune cells.Here, we model cancer-immune coevolution, wherein specialised effector cells respond to the presence of randomlyaccumulated neoantigens in a growing cancer population.We uncover a variety of cancer-immune population dynamics, from the escape or extinction of the tumour to cycles characteristic of antagonistic interactions.
We find that the suppression of the tumour by the immune system depends strongly on the cancer-immune interaction parameters, as well as rates of antigenic mutation accumulation in the cancer population.Using mutational distributions, we identify selection and the strengh of cancer-immune interactions in the system, arguing for the importance of integrating population-and single cell-level data, especially in the context of informing immunotherapeutic practices with model predictions.
Instead of encapsulating the immune impact into a selection parameter, which assumes the effector population reacts immediately and perfectly to any new antigenic mutation [11,36], we model the explicit interactions between cancer and immune cells.Our model unveils effector population dynamics during burst-like responses to growing subclones of the cancer population that possess antigenic mutations with high immunogenicity (see Figure 2).These immune population spikes serve to eliminate specific mutations from the cancer population.
Via this selective killing of cancer cells, a process known as immunoselection [15], the immune system moulds the genetic landscape of the tumour in ways that are identifiable via sequencing data.For instance, the decrease in average antigenic mutations per cell (see Figure 4) and the truncation of the high-burden tail of the MBD (see Figure 5) demonstrate that cells with more mutations face stronger negative selection by immune response and are thus pruned from the population.It should be noted, however, that for certain cancers the neoantigenic landscape has been found to be only weakly impacted by cancer-immune coevolution [10].Central to this modelling challenge is understanding the mutational process itself.
A fraction of somatic mutations arising in cancer populations give rise to the presentation of neoantigens [58].
This is a random process, wherein high mutational loads do not necessarily correspond to high antigenicities [15].
This implies that careful consideration of the mutational burdens of cells as well as the antigenicities of individual mutations is crucial to understanding the resulting evolutionary dynamics of the system.The total mutational burden of the tumour, however, is not a sufficient predictor of response to treatment unless mutations that have escaped are taken into consideration [12].While the SFS and the single-cell MBD can inform and help quantify selection (see Figures 5 and A8), more work needs to be done to understand the impact of different mechanisms of immune escape on genetic data.
We demonstrate in Section A1 of the Appendix that a higher average SFS increases the growth of effector cells: indeed, the smaller the net impact on the effector population of cancer-immune interactions (that is, the average of α i − γ i ), the larger the mean of the SFS must be to support a growing effector population (see equation (A5)).Higher intratumour heterogeneity, captured by a higher mean of the SFS, not only translates into a more effective immune response, but also to better outcomes when treated with immunotherapy [36].
Informing treatment is one of the principal tenets of mathematical modelling in oncology [32,59].Some immunotherapies increase the ability of effector cells to kill cancer cells [15], while others, termed immune checkpoint inhibitors, reactivate immune predation in the case of antigenic mutations having escaped detection [12].The neoantigens accumulated after escape thus work against the cell once immunotherapy renders the cell visible to the immune system once more, though immune evasion may still impede immunotherapy [33].
If, however, the antigenicities of the cancer cells are low due to immunoselection, the tumour will be less likely to respond well to immune checkpoint inhibition [12].
Few models have explored the relative advantages of different changes in tumour-immune interactions, which represent the impact of immunotherapies discussed above [45].Wilkie & Hahnfeldt, for instance, demonstrated that resistance to immune predation plays a smaller role than effector recruitment [16].Our results show that these relative advantages are highly dependent on the system itself: in Figure 3D, moving upwards from point c (decreasing the resistance to predation) has little impact on outcome, whereas moving to the right (increasing recruitment) changes the outcome drastically.Conversely, at point e, we notice the opposite effect: a change in predation resistance impacts the outcome but a change in recruitment does not.Importantly, only treatments that transform system parameters can succeed in a robust fashion, since only changing the state will still result in the same equilibria as before [59].
Limitations exist when trying to model the cancer-immune system.When the (antigenic) mutation rate is low, the fraction of the tumour visible to the immune system is too (see Figure A6).The cancer population dynamics are therefore largely neutral [11], though our model reveals complex effector dynamics.One can also assume a certain antigenic proportion in the tumour before immune recognition [36], or address the growth-threshold conjecture, which states that the immune system will respond to a large enough tumour growth rate, rather than a certain tumour size [33,60].The situation can also be further complicated by explicitly considering the composite state of an effector cell in the process of killing its cancer prey as a new conjugate type in the model, as has recently been done by Yang et al., who demonstrated its possible impact on the resulting dynamics, including on the outcome and its time scale [61].
By employing an individual-based model, we can compare expected mutational distributions with corresponding genetic data.The signatures of selection and strength of cancer-immune interactions in the system along with a mechanistic knowledge of these interactions may help inform us of a tumour's evolutionary history, along with its immunotherapeutic potential.
where we have written ∆ = ⟨α i − γ i ⟩ and identified that ⟨C i ⟩ is the average number of cells with a given mutation, which is the same as the average appearances of a mutation, or U/M .Thus, (A3) becomes Ignore passive recruitment (that is, the M B term in (A4)) for the moment; we will later call this the active recruitment regime, where interactions with target cancer cells contribute more to the effector population than the drift term B. A conservative condition for the effector population to be growing, dE/dt > 0, is here This inequality can be interpreted in the following way: if the net impact of effector-cancer interactions ∆ (assumed to be positive, as otherwise the effector population always decreases, barring passive recruitment effects which we assume to be small) is small, the SFS must be heavy-tailed (that is, more heterogeneous) for the effector population size to grow.

A2 The first mutant clone
Suppose an antigenic mutation arises with antigenicity A 1 and immunogenicity I 1 in a single progenitor cell within a population of neutral cancer cells.We want to inspect the ensuing clone with population C 1 and the corresponding effector population E 1 .Write b and d for the per-capita birth and (intrinsic) death rates (where we incorporate the immune escape rate into d for simplicity) of the cancer cells, regardless of mutational load, and the net growth rate (which we will always consider to be positive) r = b − d > 0. Now, the clone's population satisfies where we have written β 1 = β 0 A 1 for the killing rate of the cancer-effector interactions.From (A1), the effector population E 1 satisfies Label the effector population where the dominant recruitment mode transitions from passive to active as E 1 = E * .If we consider the coarse approximation of E * being a strict transition (that is, where all recruitment when E 1 < E * is passive and all recruitment when E 1 > E * is active; note that we will almost always be looking at populations far from E * , so we will not need to consider populations E 1 ≈ E * ), then we are interested in two time regimes: t ≪ E * /B 1 (the time at which we expect E * cells to have been recruited passively) and If the mutation arises at time t = 0, it takes 1/B 1 for the first effector cell to appear.At this point, in the deterministic framework, there are C 1 (t = 1/B 1 ) = e r/B1 cancer cells.

A3 Passive recruitment regime
For t ≪ E * /B 1 we have E 1 = B 1 t, which allows us to solve (A6) for C 1 (t): We also know that the cancer cell population at the transition time will be bounded from above: since the right-hand side of (A9) is the population supposing the effector cells have not killed any of the cancer cells.
Recall from (A7) that the transition between recruitment types occurs when which gives an equation for E * that can be numerically (though not analytically) solved: A4 Active recruitment regime When t ≫ E * /B 1 , we ignore the B 1 term of (A7) and are thus left with what appear to be predator-prey equations for E 1 and C 1 , as long as ∆ 1 is constant in C 1 (which we will see is not necessarily the case).
Because the active recruitment α 1 is a type-II functional response rather than a type-I bilinear response (like the inhibition/exhaustion term γ 1 ), as shown in Figure A1, we should subdivide this regime further into the intervals delineated by the crossover points . The net difference ∆ 1 (solid pink line) between active recruitment α 1 (dotted blue line) and inhibition/exhaustion γ 1 (dashed red line) interactions between effector and cancer cells splits the cancer cell population C 1 into three natural regions: [0, C * ), [C * , C * * ) and [C * * , ∞).Note that the active recruitment α 1 is a type-II functional response (reaching a maximum of h −1 α as C 1 → ∞), whose shape more resembles the dotted pale blue line than our piecewise linear approximation shown with a thicker line; this also translates to the thinner solid pink line for When C 1 > C * * , ∆ 1 < 0, so the effector population will decrease no matter what, barring passive recruitment.This means that if a subclone of cancer cells with a certain antigenic mutation can grow to a specific density, then it will successfully evade this immune response, since its corresponding effector population will not be able to grow to keep up with it.This, however, is restricted only to the effector type corresponding to the mutation that originated the clone; other antigenic mutations will arise in subclones, and therefore the clone overall will not have completely evaded the immune system.When C 1 < C * (especially for C 1 ≪ C * , when α 1 is well approximated by a linear function), we have a constant of motion from the predator-prey system: (Though since the passive recruitment contributes a drift term, V is not identically constant.)From standard linear stability analysis we also have equilibrium points at (C 1 , E 1 ) = (0, 0) and (D 1 /∆ 1 , r/β 1 ).Beyond the drift term that increases the population of the effectors (notably at the latter equilibrium point of the simplified system), the aforementioned effect of new antigenic mutations within the C 1 clone giving rise to other predating effector populations is also obscured by our simplifications.These both contribute to decrease the population The opposite effect is seem in the interval C * < C 1 < C * * , when ∆ 1 is positive but decreasing: the effector population E 1 may still grow, but a positive contribution to the equilibrium population of cancer cells exists in the form of the decreasing effectiveness of the effector cells.

A5 Expected mutational distributions
We begin by stating known expressions for the expected SFS and MBD in the case of neutral evolution, which we used for plots in Figures 4 and 5. From [57], the expected neutral SFS at large population size K with mutations accumulating at rate µ, conditional on survival (which we do by filtering out the realisations that do not make it to population K) is given by We take µ = λ(1 − p a ) for neutral mutations, and µ = λp a for antigenic mutations.
Under the same assumptions, the expected neutral MBD can be calculated as follows [49]: where the sum over ℓ spans at most the number of events (births or deaths), though the sum converges faster.We have used a continuous-time approximate version from [63] for the expected division distribution D ℓ instead of the discrete-time expression from [49] for ease of computation, along with the approximation that the population grows to K in log K/(b − d) time, which arises by inspection of the deterministic case of exponential growth.

A6 Neoantigen burden distribution
The immune system reacts to the antigenicities of the antigenic mutations possessed by a cancer cell rather than simply its (antigenic) mutational burden.To capture this nuance, we define the neoantigen burden N ℓ as the number of cells with cumulative antigenicity i A i (where the sum spans the antigenic mutations i within a cell) in [ℓ, ℓ + 1).Taken together, the elements N ℓ form the neoantigen burden distribution (NBD).
Following the derivation of the MBD from the division distribution in [49], we write A i,j,k ∼ Exp(1) for the antigenicity of the ith antigenic mutation of the jth cell (for some labelling of cells 1 ≤ j ≤ B k ) having a mutational burden of k and 1 A for the indicator function that is 1 on the set A and 0 otherwise.Conditioning on knowledge of the antigenic MBD {B k } M k=1 , we find A sum of k exponentially-distributed random variables with means ν obeys an Erlang distribution, which has probability density function We take expectations of both sides of equation (A14) (where ν = 1) and use the law of total expectation ] and equation (A15) to find   A2A-B) and when λ = 10 (Figures A2C-D).Both low and high immune impact parameter sets are plotted, to demonstrate that much like the MBD (see Figures 5 and A8), the tail of the distribution is shortened in the case of higher selection.This is because cells with higher antigenicities (likewise for those with higher immunogenicities) face stronger negative selection from the effector population. .Neoantigen burden distribution (NBD) (yellow lines) and the corresponding distribution for immunogenicity (purple lines) averaged over 100 realisations, when λ = 1 (A-B) and λ = 10 (C-D), for both low (dotted lines) and high (dashed lines) immune effectiveness.All parameter values not specified here are listed in Table 1 of the main text.

A7 Early stochastic dynamics
Similar to the deterministic discussion of Sections A3 and A4, passive and active recruitment regimes can also be identified in the stochastic model of the main text.On average, an effector population of type i = 1 will have passive capacity B/D: that is, at a population of E 1 = B/D, its passive recruitment rate and death rate are equal, and thus ignoring its interactions with its target population C 1 , E 1 will hover around this value.Once C 1 grows to D/(α 1 − γ 1 ), however, the active recruitment rate will equal the death rate.We can thus qualitatively describe the early dynamics of a cancer clone as follows: the mutation arises; the effector population grows to B/D; the clone grows to D/(α 1 − γ 1 ), at which point the effector population increases to properly (that is, via active recruitment) combat the threat.A4E-J).In the upper row, a blue histogram depicts the end times for the realisations, along with the average time T K to reach population size K = 3 × 10 4 , when applicable.(As seen in Figure A4I, all realisations go extinct, so no such T K is defined.)The red lines in the second row refer to the number of allowed antigenic mutations possessed by a cancer cell to consider it neutral.That is, the line labelled by 3 means that cells with 3 antigenic mutations are deemed neutral, whereas those with 4 antigenic mutations are considered antigenic and thus contribute to the proportion of the tumour that is antigenic.Note that since the immunogenicities and antigenicities are drawn from exponential distributions with mode 0, it is likely that many purportedly antigenic mutations i have very weak immunogenicity I i ≈ 0 and antigenicity A i ≈ 0. In this case, they would interact very rarely with the immune system (being effectively neutral), while still contributing to a count of cancer cells carrying antigenic mutations.
In Figure A4, the stopping time for the realisations T end was made large so to establish a good estimate on the average time T K that the cancer cell population (in the non-suppressed case) reaches the stopping population size K. Due to the stochasticity of the model, there is a distribution over end times for each realisation even for a single set of parameters.It is clear from the difference in values of the mean T K for Figures A4A and A4C (likewise for Figures A4E and A4G) that the immune system plays a role in slowing the cancer growth, since avg T A K < avg T C K (and avg T E K < avg T G K ). the second (third) row allows one (two) antigenic mutation(s) before a cell is considered antigenic.

A12 Genetics
When hyper-mutated tumours are considered (λ = 10), the genetic footprints of selection are more pronounced.
Figure A7 depicts the average number of antigenic mutations per cancer cell and Figure A8 shows the SFS and MBD in this regime.for all non-extinct realisations, with a mean given by a dashed vertical black line.In the second row, red lines with increasing paleness, labelled by i = 0, 1, 3 or 5, depict the proportion of cancer cells contain i or more antigenic mutations.All parameter values not specified here are listed in Table 1 of the main text.

C
Effector cell base dynamics.

E
Interactions between cancer and effector cells.

Figure 1 .
Figure 1.Cancer-immune stochastic model.A Cancer cells (red) stochastically divide and die with rates b andd, respectively (here ∅ represents no cell).B During a division event, daughter cells inherit all of their mother cell's antigenic and neutral mutations, depicted by numbers (where underlined numbers are antigenic).Cells carrying antigenic mutations have probability p e to escape and become neutral, as shown in the lower daughter cell with mutation 2. Each daughter cell also acquires a random number (drawn from a Poisson distribution with mean λ) of new mutations, where each mutation is antigenic with probability p a and neutral with probability 1 − p a .C For each antigenic mutation i present in the system, a corresponding effector cell population E i exists (blue), which grows with constant rate B and shrinks with per-capita death rate D. D Antigenic mutations in cancer cells (such as i and k) display unique neoantigens at the cell surface, whereas neutral mutations (such as j) do not.The neoantigens can be identified by specialised effector cells, which can only interact with the corresponding cancer cells.E When a cancer cell carrying the antigenic mutation i meets an effector cell of type i, three outcomes are possible: active recruitment of another effector cell of type i with rate α i , killing of the cancer cell with rate β i and inhibition/exhaustion of the effector cell with rate γ i .The expressions for the rates are found in equation(1).

10 Figure 3 .
Figure 3. Outcome heat maps for λ = 1 (A-C) and λ = 10 (D-F) for tumours interacting with an immune system characterised by cancer cell-effector cell interaction parameters α 0 (active recruitment of effectors) and β 0 (killing of cancer cells).A & D Proportion of suppressed tumours.B & E Proportion of extinct tumours.C & F Proportion of slow-growing tumours: tumours that are suppressed but do not go extinct.Points a and b (c and d for λ = 10) label parameter sets of low and high immune effectiveness, respectively.The red dashed line denotes γ 0 (here γ 0 = 10 −3 ), the rate of inhibition/exhaustion (see Section A4 for further discussion); see the Discussion section for details of point e.All parameter values not specified here are listed in Table1.

Figure 4 .
Figure 4. Average number of antigenic mutations per cell (solid oranges lines) for several representative realisations when λ = 1.Theoretical prediction for the accumulation of neutral mutations per cell in an exponentially-growing population shown in grey dashed line.A Low immune effect: α 0 = 0.002 and β 0 = 0.001.B High immune effect: α 0 = 0.03 and β 0 = 0.3.(Parameter sets chosen as points a and b from Figure 3A.)

Figure 5 .
Figure 5. Genetic markers of selection and cancer-immune interactions for λ = 1.We present the SFS of neutral (A & E) and antigenic mutations (B & F) and the MBD of neutral (C & G) and antigenic mutations (D & H).Left-column panels depict low immune effectiveness (α 0 = 0.002 and β 0 = 0.001) and right-column panels depict high immune effectiveness (α 0 = 0.03 and β 0 = 0.3).Black dashed lines are the theoretical predictions in the absence of an immune response.The dashed vertical lines represent the means of the MBDs in panels C, D, G & H. Results are averaged over 100 realisations, and all parameter values not specified here are listed in Table1.

)
This allows for a flexible definition of the NBD: we have chosen to discretise at integer intervals [ℓ, ℓ + 1), but could have kept it continuous with probability density function k E [B k ] t k−1 e −t /(k − 1)!.

Figure
FigureA2depicts the average NBD (along with the corresponding construction for the immunogenicity) over 100 realisations, when λ = 1 (FiguresA2A-B) and when λ = 10 (FiguresA2C-D).Both low and high immune Figure A2.Neoantigen burden distribution (NBD) (yellow lines) and the corresponding distribution for immunogenicity (purple lines) averaged over 100 realisations, when λ = 1 (A-B) and λ = 10 (C-D), for both low (dotted lines) and high (dashed lines) immune effectiveness.All parameter values not specified here are listed in Table1of the main text.

For most parameter choicesFigure A3 .
Figure A3.Genetic evidence of selection for a faster-acting immune system.A-B Average number of antigenic mutations per cell in solid oranges lines for several representative realisations when λ = 1.Theoretical prediction for the accumulation of neutral mutations per cell in an exponentially-growing population shown in grey dashed line.C-D Simulated SFS for antigenic mutations averaged over 100 realisations when λ = 10 (orange points), along with the theoretical predictions (black dashed lines) in the absence of an immune response.A & C Low immune effect: α 0 = 0.002 and β 0 = 0.001.B & C High immune effect: α 0 = 0.03 and β 0 = 0.3.(Parameter sets chosen as points a and b from Figure 3A.)

Figures
Figures A5A and A5B depict the extinction times for simulations run with λ = 1 and λ = 10, respectively.Figures A6A-C and A6D-F depicts what proportion of tumour cells carry antigenic mutations for simulations run with λ = 1 and λ = 10, respectively.The first row allows no antigenic mutations in neutral cancer cells;

10 Figure A6 .Figure A7 .Figure A8 .
Figure A6.Heat maps depicting the tumour composition (that is, the proportion of tumour cells that are antigenic) for λ = 1 (A-C) and λ = 10 (D-F).The first row is where one antigenic mutation held by a cell makes the cell antigenic; the second and third rows allow for one and two antigenic mutations (respectively) to be possessed by a cancer cell while still considering it neutral.Points a and c (respectively b and d) label parameter sets of low (respectively high) immune effectiveness.All parameter values not specified here are listed in Table of the main text.
Probability of escape of an antigenic cancer cell into neutrality C iPopulation of cancer cells carrying antigenic mutation i, with total cancer cell population C E i simulation time).Because the model is stochastic, different outcomes (suppression or not) may arise from a same set of parameters and initial conditions.

Table 1
of the main text.