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 evolution-ary 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, formulated by Burnet as the immunosurveillance hypothesis, 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 positive feedback loop 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 deathligand 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 size [15, 23]; immune checkpoint blockade therapies restore the effector response to immune-escaped 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 A [41], as shown in Figure 1B. Similar to Lakatos et al., we focus on the exomic region [11], where the value of A 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]. Here, we write the antigenic mutations of a cancer cell (labelled by𝓁) as M_a,𝓁 and the neutral mutations as M_n,𝓁, where M_𝓁 = M_a,𝓁 ⋃ M_n,𝓁 is the set of all mutations carried by a cancer cell. 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), rendering their antigenic mutations sets empty. It is important to note that immune escape is not necessarily a permanent state; cells that have undergone escape may subsequently acquire new antigenic mutations in future divisions, thereby regaining susceptibility to immune detection.

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, as in gene-for-gene coevolution [37]. We will write (E_i) to denote an effector cell of type i and E_i to represent the corresponding population size (and, in an abuse of notation, occasionally the population itself, when this will not cause too much confusion). Whenever mutation i exists within the population, effector cells of type i are passively recruited from the body at constant rate , for Ø the absence of a cell, 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 effector cells of type i are recruited 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. We will write (C_i) to denote a cancer cell possessing an antigenic mutation i and C_i for the corresponding population size. Interactions between cancer cells (C_i) and effector cells (E_i) have three possible outcomes: active recruitment of another effector cell (E_i) with rate α_i, killing of the cancer cell (C_i) with rate β_i and inhibition/exhaustion of the effector cell (E_i) with rate γ_i (see Figure 1E). The precise description of the system can be described by the following set of microscopic reactions:

where , and α₀, β₀, γ₀, h_α are constants for the entire tumour. For simplicity, we will call α₀, β₀ and γ₀ the recruitment, killing and inhibition rates, respectively, even though they are not strictly rates. We have denoted the cell carrying a set of mutations M_𝓁, which can be partitioned into antigenic mutations M_a,𝓁 and neutral mutations M_n,𝓁, by . Upon division, each daughter cell inherits the full set of mutations from the mother cell and independently acquires additional new mutations, denoted by M ⁰ and M ⁰⁰ for the two daughter cells, respectively. The numbers of newly acquired antigenic and neutral mutations for each daughter cell are drawn independently from Poisson distributions with means p_aλ and (1 - p_a)λ, respectively. With probability 1 - p_e, the daughter cell retains the antigenic and neutral mutations of the mother cell. With probability p_e, the daughter cell undergoes immune escape, and all mutations (including those inherited and newly acquired) are considered neutral, resulting in empty antigenic mutation set.

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 population sizes 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; this is because there is an upper bound to how quickly new effector cells can be recruited [47, 48].

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_a denotes the total number of antigenic mutations across all cancer cells [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 size exceeds a threshold K, which we call no suppression; goes extinct ; or the cancer cell remains below K at time T_end, which we call slow-growing. We call the latter two outcomes—extinction and slow-growth of the tumour cell population—suppression. While our system is stochastic, different outcomes may arise among individual realisations under the same parameter set and initial condition. Thus, we measure the proportion of simulations ending in a given outcome under different parameter values.

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 recruitment and killing rates α₀ and β₀, 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, following classic in-phase cycles with the cancer population as seen in other antagonistic systems [50].

Figure 2.

Higher mutation rates lead to the appearance of more mutants, thus potentially more new types of cancer– immune interactions. Based on previous predator–prey studies [51, 52], we developed a method to quantify the number of cancer–immune cycles at low and high mutation rates. Different from typical predator–prey systems where the two antagonistic species fluctuate between relatively stable ranges of population sizes [51], the abundances of effector and cancer cells in our system often have increasing trends, as seen in Figures 2A and 2C. Consequently, the phase portraits of the abundances of cancer effectors show stochastic and irregular behaviour (see Figure A4), rather than having a closed or open oval shape as in predator–prey systems. This makes it hard to quantify whether a cycle is in-phase or out-of-phase using the shape of the phase portraits as in predator–prey systems [53, 54], although visually out-of-phase cycles are rare in our simulations. Instead, we develop a method to quantify the number of cancer–immune cycles in our simulations by tracking the directional changes in phase portraits, validated by using a non-evolving stochastic predator–prey system as a control (see Figure A5). As expected, the number of cycles increases when the mutation rate is higher (see Figure A6). The majority are counter-clockwise cycles, where the cancer population increases first followed by the increase of the effector population (see Figure A7). However, we also observe a small fraction of clockwise cycles, especially when mutation rate is higher. Clockwise cycles have been observed in various predator–prey systems, and could arise as a consequence of coevolution [55].

Because of the specialised nature of our model, these cyclic dynamics arise when a single antigenic mutation i causes the rapid active recruitment of the corresponding effector cells (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 A4 of the Appendix. We validate this by inspecting the coevolution of antigenicity and immunogenicity in the corresponding single realisations. In Figures 2A & 2C, we can see that the effector spikes arise due to the emergence of one or several mutations that has a large immunogenicity I_i (Figure A3) dominating the average immunogenicity in the effector population (Figure 2B& 2D)). Unexpectedly, in addition to this elevated immunogenicity, Figures 2B& 2D as well as Figure A3 imply that these spikes also arise for mutations with low antigenicity. These fluctuations in ⟨I⟩ and ⟨A⟩ are coupled to the spikes in the population size (Figures 2 and A3); 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 more cycles of spikes not only in the cancer–effector dynamics but also 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 Figure A10A–C versus Figure A10D–F of the Appendix. It also worth noting that in this case most suppressed tumours go extinct at early times (Figure A9A) 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. 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 recruitment and killing rates α₀ 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^-1) 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 and Hahnfeldt [16].

Genetic markers of selection

Stochastic mutation accumulation in an exponentially growing population has been widely studied [56, 57]. 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 [58], as plotted in Figure 4. (Note that a correction for large times has recently been shown to apply [59, 60], 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 be 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 A12 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 [61]. As discussed in Section A1 of the Appendix, the expected neutral distributions of the SFS and the MBD have been solved [49, 62], 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 (A1) and (A2) 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 A13 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 A5.

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 earlystage small cell population sizes 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 strength 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 [63]. 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 A13), more work needs to be done to understand the impact of different mechanisms of immune escape on genetic data.

Informing treatment is one of the principal tenets of mathematical modelling in oncology [32, 64]. 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 and 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 [64].

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 A10). 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, 65]. 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 [66].

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.

Additional files

Supplementary text and figures