Figures and data in Cancer–immune coevolution dictated by antigenic mutation accumulation

Figures
Tables
Additional files

19 figures, 1 table and 1 additional file

Figures

Figure 1

Download asset Open asset

Cancer–immune stochastic model.

(A) Cancer cells (red) stochastically divide and die with rates $b$ and $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 a 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), whose size 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.

Figure 2

Download asset Open asset

Single stochastic realisations for mutation rate $λ = 1$ (**A, B**) and $λ = 10$ (**C, D**).

(**A, C**) Population dynamics, where red and blue lines depict total cancer and effector cell populations. (**B, D**) Average immunogenicity (yellow line) and antigenicity (orange line) in the effector population for the corresponding realisation. Vertical dashed grey lines in all panels indicate the timing of effector population spikes, which are consistent with the vertical dashed lines in the Muller plots (Appendix 1—figure 3) showing the dynamics of individual effector phenotypes over time for the same realisations. Interaction parameters: recruitment rate $α_{0} = 0.03$ and killing rate $β_{0} = 0.3$ (**A, B**) and $α_{0} = 0.005$ and $β_{0} = 0.01$ (**C, D**).

Figure 3

Download asset Open asset

Outcome heat maps for mutation rate $λ = 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 inhibition rate $γ_{0}$ (here $γ_{0} = 10^{- 3}$ ); see the Discussion section for details of point $e$ . All parameter values not specified here are listed in Table 1.

Figure 4

Download asset Open asset

Average number of antigenic mutations per cell (solid orange lines) for several representative realisations when $λ = 1$ .

Theoretical prediction for the accumulation of neutral mutations per cell in an exponentially growing population is shown in grey dashed line. (A) Low immune effect: recruitment rate $α_{0} = 0.002$ and killing rate $β_{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

Download asset Open asset

Genetic markers of selection and cancer–immune interactions for $λ = 1$ .

Panels show the site frequency spectrum (SFS) for neutral (**A, E**) and antigenic mutations (**B, F**), and the single-cell mutational burden distribution (MBD) for neutral (**C, G**) and antigenic mutations (**D, H**). Left-column panels correspond to low immune effectiveness ( $α_{0} = 0.002$ , $β_{0} = 0.001$ , point a in Figure 3), while right-column panels correspond to high immune effectiveness (recruitment rate $α_{0} = 0.03$ , killing rate $β_{0} = 0.3$ , point b in Figure 3). Black dashed lines indicate theoretical predictions under neutral evolution (no immune response). Solid and dashed vertical lines in panels C, D, G, and H denote the means of the simulated and theoretical MBDs, respectively. The $W_{1}$ value represents the Wasserstein distance between the simulated and theoretical MBDs with the number of mutations $k$ is rescaled to the interval [0,1], and the number of cells $B_{k}$ is normalised to form a probability density. Results are averaged over 100 realisations; all other parameter values are given in Table 1.

Appendix 1—figure 1

Download asset Open asset

Genetic evidence of selection for a faster-acting immune system.

(**A, B**) Average number of antigenic mutations per cell in solid orange 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).

Appendix 1—figure 2

Download asset Open asset

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 (solid lines) immune effectiveness.

All parameter values not specified here are listed in Table 1 of the main text.

Appendix 1—figure 3

Download asset Open asset

Muller plots for two representative realisations from Figure 2.

(A, B) correspond to Figure 2A, B and C, D, respectively. Vertical dashed grey lines in all panels refer to the timing of spikes in the effector population in Figure 2, also indicated by vertical dashed grey lines. The colour indicates the antigenicity and immunogenicity of the antigen corresponding to each effector type, as shown in the colour maps on the right. White regions indicate the absence of effector cells, while grey represents all rare effector types whose abundance remains below 5%. Interaction parameters: $α_{0} = 0.03$ , $β_{0} = 0.3$ (A); $α_{0} = 0.005$ , $β_{0} = 0.01$ (B).

Appendix 1—figure 4

Download asset Open asset

Phase trajectories of two realisations in Figure 2.

(**A, B**) are corresponding to Figure 2A, B and C, D, respectively. Interaction parameters: $α_{0} = 0.03$ , $β_{0} = 0.3$ (A); $α_{0} = 0.005$ , $β_{0} = 0.01$ (B).

Appendix 1—figure 5

Download asset Open asset

Validation of our cycle counting method in a non-evolving prey–predator system.

(A) Population dynamics of the deterministic system (dark lines) with 20 stochastic realisations (pale lines). (B) Inferred number of predator–prey cycles using the same methodology as in Appendix 1—figure 6. Model parameters, from Equation 9: prey birth rate = 1, prey death rate = 0.1, predator death rate = 0.1, predation rate = 0.001, and predator reproduction rate = 0.0001. Initial population sizes of the prey and the predator were 1000 and 700, respectively, with a maximum simulation time of 100. The histogram is constructed from 500 independent realisations.

Appendix 1—figure 6

Download asset Open asset

Histograms depicting the inferred number of cycles in the population dynamics.

(**A, B**) Distribution of cycle counts for simulations with, respectively, low ( $λ = 1$ ) and high ( $λ = 10$ ) mutation rates. Each histogram is constructed from 200 independent realisations. The interaction parameters are $α_{0} = 0.03$ and $β_{0} = 0.3$ (A); $α_{0} = 0.005$ and $β_{0} = 0.01$ (B).

Appendix 1—figure 7

Download asset Open asset

Histograms depicting the inferred number of irregular (clockwise) cycles in the population dynamics.

(**A, B**) Distribution of cycle counts for simulations with, respectively, low ( $λ = 1$ ) and high ( $λ = 10$ ) mutation rates. Each histogram is constructed from 200 independent realisations. The interaction parameters are $α_{0} = 0.03$ and $β_{0} = 0.3$ (A); $α_{0} = 0.005$ and $β_{0} = 0.01$ (B).

Appendix 1—figure 8

Download asset Open asset

Cancer–effector population dynamics (first row) and tumour antigenicity (second row) for five different sets of parameters.

(**A, B**) $λ = 1$ , $α_{0} = 0.002$ and $β_{0} = 0.001$ . (**C, D**) $λ = 1$ , $α_{0} = 0.03$ and $β_{0} = 0.3$ . (**E, F**) $λ = 10$ , $α_{0} = 0.002$ and $β_{0} = 0.001$ . (**G, H**) $λ = 10$ , $α_{0} = 0.005$ and $β_{0} = 0.01$ . (**I, J**) $λ = 10$ , $α_{0} = 0.03$ and $β_{0} = 0.3$ . In the first row, red lines depict cancer cell population sizes for several representative realisations. Blue histograms above show the final time $T_{K}$ where the population size reached $K = 3 \times 10^{4}$ 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.

Appendix 1—figure 9

Download asset Open asset

Extinction time heat maps for $λ = 1$ (A) and $λ = 10$ (B). All parameter values not specified here are listed in Table 1 of the main text.

Appendix 1—figure 10

Download asset Open asset

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 (A & D) is where one antigenic mutation held by a cell makes the cell antigenic; the second (B & E) and third rows (C & F) allow for one and two antigenic mutations (respectively) to be possessed by a cancer cell while still considering it neutral. Points 𝑎 and 𝑐 (respectively, 𝑏 and 𝑑) label parameter sets of low (respectively, high) immune effectiveness. All parameter values not specified here are listed in Table 1 of the main text.

Appendix 1—figure 11

Download asset Open asset

Genetic markers of selection: SFS (**A, B**) and MBD (**C, D**) for no immune effect case (all interaction parameters set to 0) averaged over 100 realisations when $λ = 1$ , along with the theoretical predictions (black dashed lines) in the absence of an immune response.

Green data represents neutral mutations (**A, C**) and orange data represents antigenic mutations (**B, D**), with vertical lines representing the means of the distributions for MBDs in panels (C, D). All parameter values not specified here are listed in Table 1 of the main text.

Appendix 1—figure 12

Download asset Open asset

Average number of antigenic mutations in solid orange lines for several representative realisations when $λ = 10$ .

Theoretical prediction for the accumulation of neutral mutations in an exponentially growing population is shown in grey dashed line. (A) Low immune effect: $α_{0} = 0.002$ and $β_{0} = 0.001$ . (B) Middling immune effect: $α_{0} = 0.005$ and $β_{0} = 0.01$ . (C) High immune effect: $α_{0} = 0.03$ and $β_{0} = 0.3$ . (Parameter sets for (A) and (C) chosen as points $c$ and $d$ from Figure 3C, respectively.)

Appendix 1—figure 13

Download asset Open asset

Genetic markers of selection: SFS (**A, B**, **E, F**) and MBD (**C, D**, **G, H**) for low (left column; $α_{0} = 0.002$ and $β_{0} = 0.001$ , point c in Figure 3) and high (right column; $α_{0} = 0.005$ and $β_{0} = 0.01$ , point d in Figure 3) immune effectiveness averaged over 100 realisations when $λ = 10$ , along with the theoretical predictions (black dashed lines) in the absence of an immune response.

Green data represents neutral mutations (A, C, **E, G**) and orange data represents antigenic mutations (B, D, **F, H**), with vertical lines representing the means of the distributions for MBDs in panels (C, D, **G, H**). All parameter values not specified here are listed in Table 1 of the main text.

Appendix 1—figure 14

Download asset Open asset

Heat maps illustrate the differences between the simulated and theoretical single-cell mutation burden distributions when $λ = 1$ .

The colour gradient represents the Wasserstein distance, with the number of mutations $k$ rescaled to the interval [0,1]. The left panel displays how the difference between simulated and theoretical predicted neutral MBDs varies with interaction parameters, while the right panel shows the same for antigenic mutations. Point a and b are the same in Figure 3 and all other parameter values can be found in Table 1 in the main text.

Tables

Table 1

Notation and baseline parameter values.

Symbol	Baseline	Description
$b$ , $d$	1 , 0.1	Birth and death rates of a cancer cell
$λ$	1 , 10	Mean number of exomic mutations acquired per division per daughter cell
$p_{a}$	0.075	Probability of a mutation being antigenic (rather than neutral)
$p_{e}$	10^–6	Probability of escape of an antigenic cancer cell into neutrality
$C_{i}$		Population of cancer cells carrying antigenic mutation $i$ , with total cancer cell population size $C$
$E_{i}$		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 immunogenicity, 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_{a, ℓ}, M_{n, ℓ}$		Set of antigenic and neutral mutations, respectively, carried by a cancer cell $ℓ$
$M_{a}, M_{n}$		Total number of all antigenic and neutral mutations: $M_{a} = \| ⋃_{ℓ} M_{a, ℓ} \|$ and $M_{n} = \| ⋃_{ℓ} M_{n, ℓ} \|$
$U$		Total antigenic mutational occurrences: $U = \sum_{j = 1}^{C} j S_{j} = \sum_{k = 1}^{M_{a}} k B_{k}$
$K$	3 × 10⁴	Ending cancer population size for stochastic simulations
$T_{end}$	16	Maximum time for a simulation

Additional files

MDAR checklist: https://cdn.elifesciences.org/articles/103970/elife-103970-mdarchecklist1-v1.pdf
Download elife-103970-mdarchecklist1-v1.pdf

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Article PDF

Open citations (links to open the citations from this article in various online reference manager services)

Mendeley

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Long Wang
Christo Morison
Weini Huang

(2025)

Cancer–immune coevolution dictated by antigenic mutation accumulation

eLife 14:RP103970.

https://doi.org/10.7554/eLife.103970.3

Figures

Cancer–immune stochastic model.

Single stochastic realisations for mutation rate $λ = 1$ (A, B) and $λ = 10$ (C, D).

Outcome heat maps for mutation rate $λ = 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).

Average number of antigenic mutations per cell (solid orange lines) for several representative realisations when $λ = 1$ .

Genetic markers of selection and cancer–immune interactions for $λ = 1$ .

Genetic evidence of selection for a faster-acting immune system.

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 (solid lines) immune effectiveness.

Muller plots for two representative realisations from Figure 2.

Phase trajectories of two realisations in Figure 2.

Validation of our cycle counting method in a non-evolving prey–predator system.

Histograms depicting the inferred number of cycles in the population dynamics.

Histograms depicting the inferred number of irregular (clockwise) cycles in the population dynamics.

Cancer–effector population dynamics (first row) and tumour antigenicity (second row) for five different sets of parameters.

Extinction time heat maps for $λ = 1$ (A) and $λ = 10$ (B). All parameter values not specified here are listed in Table 1 of the main text.

Heat maps depicting the tumour composition (that is, the proportion of tumour cells that are antigenic) for $λ = 1$ (A—C) and $λ = 10$ (D—F).

Genetic markers of selection: SFS (A, B) and MBD (C, D) for no immune effect case (all interaction parameters set to 0) averaged over 100 realisations when $λ = 1$ , along with the theoretical predictions (black dashed lines) in the absence of an immune response.

Average number of antigenic mutations in solid orange lines for several representative realisations when $λ = 10$ .

Heat maps illustrate the differences between the simulated and theoretical single-cell mutation burden distributions when $λ = 1$ .

Tables

Notation and baseline parameter values.

Additional files

MDAR checklist

Download links

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Be the first to read new articles from eLife

Share this article

Cite this article

Cancer–immune stochastic model.

Single stochastic realisations for mutation rate λ=1 (A, B) and λ=10 (C, D).

Outcome heat maps for mutation rate λ=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).

Average number of antigenic mutations per cell (solid orange lines) for several representative realisations when λ=1.

Genetic markers of selection and cancer–immune interactions for λ=1.

Genetic evidence of selection for a faster-acting immune system.

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 (solid lines) immune effectiveness.

Muller plots for two representative realisations from Figure 2.

Phase trajectories of two realisations in Figure 2.

Validation of our cycle counting method in a non-evolving prey–predator system.

Histograms depicting the inferred number of cycles in the population dynamics.

Histograms depicting the inferred number of irregular (clockwise) cycles in the population dynamics.

Cancer–effector population dynamics (first row) and tumour antigenicity (second row) for five different sets of parameters.

Extinction time heat maps for λ=1 (A) and λ=10 (B). All parameter values not specified here are listed in Table 1 of the main text.

Heat maps depicting the tumour composition (that is, the proportion of tumour cells that are antigenic) for λ=1 (A—C) and λ=10 (D—F).

Genetic markers of selection: SFS (A, B) and MBD (C, D) for no immune effect case (all interaction parameters set to 0) averaged over 100 realisations when λ=1, along with the theoretical predictions (black dashed lines) in the absence of an immune response.

Average number of antigenic mutations in solid orange lines for several representative realisations when λ=10.

Heat maps illustrate the differences between the simulated and theoretical single-cell mutation burden distributions when λ=1.

Notation and baseline parameter values.

MDAR checklist

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Single stochastic realisations for mutation rate $λ = 1$ (A, B) and $λ = 10$ (C, D).

Outcome heat maps for mutation rate $λ = 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).

Average number of antigenic mutations per cell (solid orange lines) for several representative realisations when $λ = 1$ .

Genetic markers of selection and cancer–immune interactions for $λ = 1$ .

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 (solid lines) immune effectiveness.

Extinction time heat maps for $λ = 1$ (A) and $λ = 10$ (B). All parameter values not specified here are listed in Table 1 of the main text.

Heat maps depicting the tumour composition (that is, the proportion of tumour cells that are antigenic) for $λ = 1$ (A—C) and $λ = 10$ (D—F).

Genetic markers of selection: SFS (A, B) and MBD (C, D) for no immune effect case (all interaction parameters set to 0) averaged over 100 realisations when $λ = 1$ , along with the theoretical predictions (black dashed lines) in the absence of an immune response.

Average number of antigenic mutations in solid orange lines for several representative realisations when $λ = 10$ .

Heat maps illustrate the differences between the simulated and theoretical single-cell mutation burden distributions when $λ = 1$ .