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 immunosurveil-lance 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–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 tumour-immune 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.

Model

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).

Figure 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):

where α₀, β₀, γ₀ 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 α₀ > γ₀, 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 single-cell 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 simulation time). Because the model is stochastic, different outcomes (suppression or not) may arise from a same set of parameters and initial conditions.

Results

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 α₀ and β₀ (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 neutral-like 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.

Figure 2.

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 α₀ and the cancer killing rate β₀. Due to the stochastic nature of the model, the outcome for a given parameter set is probabilistic. As expected, for higher values of α₀ and β₀, 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.

Figure 3.

When the active effector recruitment rate is smaller than the effector inhibiting rate α₀ < γ₀ (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 α₀ 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 α₀ and β₀ (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.

Table 1.

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 β₀ ≈ 10⁻¹) above which tumours are suppressed, no matter the active recruitment rate α₀. The rest of the domain of Figure 3D then exhibits a much stronger dependence on the active recruitment rate α₀, in line with the results of Wilkie & Hahnfeldt [16].

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 of the Appendix.

Figure 4.

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.

Figure 5.

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.

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.

Code Availability

Code of implementation of our model is available at https://github.com/Bio421/CancerImmuneCoevo under the MIT license.

Acknowledgements

L.W. was supported by the International Program for Candidates, Sun Yat-Sen University. C.M. was supported by the European Union’s Horizon 2020 research and innovation programme under the Marie SkŁodowska-Curie EvoGamesPlus grant number 955708.

Appendix

This appendix is split into three sections: deterministic analysis, stochastic analysis and extended figures. We first state the ODE limit of our stochastic model and draw some conclusions about intratumour heterogeneity and its impact on effector cell population growth (Section A1). We then inspect the first mutant clone to arise in the cancer cell population (Section A2), which suggests investigating the two regimes of effector cell population growth when passive and active recruitment dominate, respectively (Sections A3 and A4). We return to the stochastic model to state theoretical distributions from the literature (Section A5). Second, we define the neoantigen burden distribution and derive its expected distribution, as well as plot the correlation between the mean of the NBD and the mean of the SFS (Section A6). We then discuss early dynamics (Section A7). Next, we investigate the effector cell population spikes and derive approximations for their expected frequency and amplitudes (Section A8). Finally, we discuss the implications of a faster immune response due to an existing effector pool in the body on the genetic evidence for selection A9. Additional figures of stochastic simulations are found in Sections A10-A12.

Deterministic analysis

A1 ODE model

Associate to the antigenic mutation labelled i an effector population E_i with passive recruitment rate B and per-capita death rate D, along with α_i − γ_i the factor denoting the net difference between active recruitment and inhibition/exhaustion interactions between effector cells of type i and cancer cells possessing mutation i, whose population we denote C_i. The ODE describing the population dynamics of E_i is thus

We will write E = ∑_i E_i for the total effector population. We will ignore neutral mutations in the ODE model, so unless otherwise specified, mutations are now assumed to be antigenic.

Because we are considering the emergence of cancer (that is, our populations are small), we model tumour growth exponentially. It is expected that the immune system is able to handle small tumours [5], so the tumour per-capita growth rate should, in principle, be bounded near populations of zero [31]. Note that the oft-used Gompertz model of tumour growth, which is unbounded at zero, was originally used for tumours whose initial sizes were on the order of thousands of cells [62].

Recall the following definitions from the main text. Let the elements of the site frequency spectrum (SFS) S_j denote the number of mutations occurring j times in the (cancer cell) population, so that M = ∑_j S_j is the total number of mutations. Let the elements of the single-cell mutational burden distribution (MBD) B_k denote the number of cells possessing k mutations, so that C = ∑_k B_k is the total cancer population. The total number of mutational occurrences (see the Supplementary Information of [49] for details) is then

so that the average appearances of a mutation (that is, the mean of the SFS) is U/M, and the average single-cell mutational burden (that is, the mean of the MBD) is U/N.

Notice that for a fixed mutation i, when α_i − γ_i is high, we expect C_i to be low, since more effector cells of type i are around (having been recruited but not inhibited nor exhausted) to deplete the population C_i. Thus we expect that cov(α_i − γ_i, C_i) ≤ 0 for all i.

Hoping to inspect the overall effector population, we sum (A1) over mutations 1 ≤ i ≤ M to find

for ⟨·⟩ the average over i. Now,

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 MB 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₁ and immunogenicity I₁ in a single progenitor cell within a population of neutral cancer cells. We want to inspect the ensuing clone with population C₁ and the corresponding effector population E₁. 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 β₁ = β₀A₁ for the killing rate of the cancer-effector interactions. From (A1), the effector population E₁ satisfies

Label the effector population where the dominant recruitment mode transitions from passive to active as E₁ = E^*. If we consider the coarse approximation of E^* being a strict transition (that is, where all recruitment when E₁ < E^* is passive and all recruitment when E₁ > 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₁ ≈ E^*), then we are interested in two time regimes: t ≪ E^*/B₁ (the time at which we expect E^* cells to have been recruited passively) and t ≫ E^*/B₁.

If the mutation arises at time t = 0, it takes 1/B₁ for the first effector cell to appear. At this point, in the deterministic framework, there are C₁(t = 1/B₁) = e^r/B1 cancer cells.

A3 Passive recruitment regime

For t ≪ E^*/B₁ we have E₁ = B₁t, which allows us to solve (A6) for C₁(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 B₁ = Δ₁E₁C₁ (where Δ₁ = α₁₋γ₁), which gives an equation for E^* that can be numerically (though not analytically) solved:

A4 Active recruitment regime

When t ≫ E^*/B₁, we ignore the B₁ term of (A7) and are thus left with what appear to be predator-prey equations for E₁ and C₁, as long as Δ₁ is constant in C₁ (which we will see is not necessarily the case). Because the active recruitment α₁ is a type-II functional response rather than a type-I bilinear response (like the inhibition/exhaustion term γ₁), as shown in Figure A1, we should subdivide this regime further into the intervals delineated by the crossover points C₁ = C^* and C₁ = C^**.

Figure A1.

When C₁ > C^**, Δ₁ < 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₁ < C^* (especially for C₁ ≪ C^*, when α₁ 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₁, E₁) = (0, 0) and (D₁/Δ₁, r/β₁). 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₁ clone giving rise to other predating effector populations is also obscured by our simplifications. These both contribute to decrease the population C₁.

The opposite effect is seem in the interval C^* < C₁ < C^**, when Δ₁ is positive but decreasing: the effector population E₁ 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.

Stochastic analysis

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 𝟙_A for the indicator function that is 1 on the set A and 0 otherwise. Conditioning on knowledge of the antigenic MBD , 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 𝔼 [X] = 𝔼 [𝔼 [X|Y]] and equation (A15) to find

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 𝔼 [B_k] t^k−1e^−t/(k − 1)!.

Figure A2 depicts the average NBD (along with the corresponding construction for the immunogenicity) over 100 realisations, when λ = 1 (Figures 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.

Figure A2.

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₁ = B/D, its passive recruitment rate and death rate are equal, and thus ignoring its interactions with its target population C₁, E₁ will hover around this value. Once C₁ grows to D/(α₁ − γ₁), 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/(α₁ − γ₁), at which point the effector population increases to properly (that is, via active recruitment) combat the threat.

For most parameter choices in the main text, B/D = O(1) and D/(α₁ − γ₁) = O(10²), so that the effector population only manages to grow to more than a few cells once its targets number around a hundred. (It is also worth mentioning that early ODE models of tumour-immune dynamics phrased in a predator-prey fashion use Kuznetsov et al.’s estimate of effector passive recruitment being O(10⁴) cells per day [44]. With O(0.1) or O(1) new antigenic mutations (and thus effector types) generated in each division, when our cancer cell population is O(10⁴), we therefore expect a passive recruitment similar to Kuznetsov et al. [44].)

Models of immune responses vary in their triggers: Chen & Xie et al. suppose a minimal proportion of antigenic cancer cells [36], the growth threshold conjecture requires a minimal cancer growth rate [64], and most ODE models as well as Lakatos et al. suppose instantaneous responses [11, 31]. While our mechanistic model does not assume frequency dependence, the previous discussion implies that the immune response implicitly requires a certain cancer density before it can grow effectively to fight.

A8 Effector cell population spikes

As identified in individual-based simulations (see Figure 2, for instance), effector subpopulations in the stochastic model can undergo rapid spiking events, before vanishing to zero when the targeted mutant (labelled by i) goes extinct (and thus the effector cells have no more prey and so are removed from the system). These spikes arise when the mutation has immunogenicity I_i much higher than the mean ⟨I⟩, as well as an antigenicity A_i lower than the mean ⟨A⟩, as otherwise the spike would not need to grow since the effector cells would be more efficient at killing. Qualitatively, we will write Î (Â) for the threshold immunogenicity (antigenicity) and above (below) which spikes may occur.

Returning to the language of predator-prey dynamics discussed in Section A4, effector spikes arise when the populations enter a cycle that ends when C_i(t) < 1 (that is, extinction of mutation i). We write i = 1 for the mutation causes the spike, since we consider it to be arising in a neutral cell. (We argue qualitatively that were other antigenic mutations carried by the cell to generate an effector spike, they already would have; thus we expect most spikes to arise in neutral or almost neutral cells.) Now, recalling the constant of motion (A11) and setting C₁ = 1, we obtain an equation for E₁ (as a function of α₁ − γ₁), the approximate (due to the stochasticity of the model as well as the simplifying assumptions above) amplitude of the spike:

which can have up to two solutions and can be numerically solved. Differentiating (A17) implicitly with respect to I₁, we obtain an ODE for E₁:

When the cancer clone is decreasing in population, the first term is negative; the second term is always negative, since α₁ is increasing in I₁ by inspection of (1). Thus, we have that dE₁/dI₁ > 0, confirming that a higher immunogenicity does increase the corresponding effector population, as expected.

Using that immunogenicities and antigenicities are drawn from exponential distributions with means 1, the rate of spikes r_spike is then

since the rate of new antigenic mutations and the probabilities ℙ(I_i > Î) and ℙ(A_i > Â) are independent.

A9 Faster immune response

Our immune system possesses an innate pool of effector cells, able to respond to threats without the timeconsuming training assumed implicitly in our specialised model of active and passive recruitment [19]. With some probability, this pool might include effector cells that can recognise newly arising neoantigens quickly. Thus, rather than having to wait for the first passively-recruited effector cell, the body could react faster to novel antigenic mutations.

To investigate the impact of this mechanism on the observed coevolution, we modified our model to include an effector cell of type i upon initiation of antigenic mutation i. Figures A3A-B show the average number of antigenic mutations per cell and Figures A3C-D show the SFS, in the case of λ = 1.

When comparing with the corresponding figures for the model used throughout the main text (Figures 4A-B, 5B and 5F, respectively), the patterns appear the same, though more evident. This reinforces that the speed of the immune response increases its effectiveness though does not change the qualitative nature of cancer-immune coevolution. However, rather than having a guaranteed effector cell that can handle a newly-arising neoantigen, there is some probability of the effector pool containing the desired specialty. Thus we expect that in the body, selection will follow the patterns shown in both cases (Figures 4 and 5 versus Figure A3), at some magnitude between the two, related to the aforementioned probability.

Figure A3.

Extended figures

A10 Population dynamics

Figure A4 shows representative population dynamics (first row), along with the corresponding composition of the tumour population (second row), when λ = 1 (Figures A4A-D) and λ = 10 (Figures 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⁴, 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 (and ).

A11 Heat maps

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; 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.

Figure A4.

Figure A5.

Figure A6.

Figure A7.

Figure A8.