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 probability pe 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 A) of new mutations, where each mutation is antigenic with probability pa and neutral with probability 1 - pa. C For each antigenic mutation i present in the system, a corresponding effector cell population Ei 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).

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 consists with the vertical dashed lines in the Muller plots (Figure A3) 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).

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.

Notation and baseline parameter values, which are used for all simulations unless otherwise specified.

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

Theoretical prediction for the accumulation of neutral mutations per cell in an exponentially-growing population shown in grey dashed line. A Low immune effect: 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.)

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), while right-column panels correspond to high immune effectiveness (recruitment rate α0 = 0.03, killing rate β0 = 0.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 W1 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 Bk is normalised to form a probability density. Results are averaged over 100 realisations; all other parameter values are given in Table 1.