Luis Almeida
email: almeida@ljll.math.upmc.fr

Chloe Audebert
email: chloe.audebert@sorbonne-universite.fr

Emma Leschiera
email: leschiera@ljll.math.upmc.fr

Tommaso Lorenzi
email: tommaso.lorenzi@polito.it

Discrete and continuum models for the coevolutionary dynamics between CD8+ cytotoxic T lymphocytes and tumour cells

Keywords: Immune competition, Coevolutionary dynamics between cytotoxic T lymphocytes and tumour cells, Patterns of phenotypic coevolution, Individual-based models, Continuum models

We present an individual-based model for the coevolutionary dynamics between CD8+ cytotoxic T lymphocytes (CTLs) and tumour cells. In this model, every cell is viewed as an individual agent whose phenotypic state is modelled by a discrete variable. For tumour cells this variable represents a parameterisation of the antigen expression profiles, while for CTLs it represents a parameterisation of the target antigens of T-cell receptors (TCRs). We formally derive the deterministic continuum limit of this individual-based model, which comprises a non-local partial differential equation for the phenotype distribution of tumour cells coupled with an integro-differential equation for the phenotype distribution of CTLs. The biologically relevant homogeneous steady-state solutions of the continuum model equations are found. The linear-stability analysis of these steady-state solutions is then carried out in order to identify possible conditions on the model parameters that may lead to different outcomes of immune competition and to the emergence of patterns of phenotypic coevolution between tumour cells and CTLs. We report on computational results of the individual-based model, and show that there is a good agreement between them and analytical and numerical results of the continuum model. These results shed light on the way in which different parameters affect the coevolutionary dynamics between tumour cells and CTLs. Moreover, they support the idea that TCR-tumour antigen binding affinity may be a good intervention target for immunotherapy and offer a theoretical basis for the development of anti-cancer therapy aiming at engineering TCRs so as to shape their affinity for cancer targets.

Introduction

Essentials of the underlying biological problem CD8+ cytotoxic T lymphocytes (CTLs) play a key role in the immune response against cancer. CTLs carry a specific receptor on their surface, the T-cell receptor (TCR), which can recognise and bind to non-self antigens expressed by tumour cells [START_REF] Messerschmidt | How cancers escape immune destruction and mechanisms of action for the new significantly active immune therapies: Helping nonimmunologists decipher recent advances[END_REF]. Each TCR recognises and binds specifically to a certain antigen (i.e. the cognate antigen) [START_REF] Coulie | Tumour antigens recognized by T lymphocytes: at the core of cancer immunotherapy[END_REF], and possibly other antigens within a certain affinity range [START_REF] Mason | A very high level of crossreactivity is an essential feature of the T-cell receptor[END_REF][START_REF] Wooldridge | A single autoimmune T cell receptor recognizes more than a million different peptides[END_REF]. This enables CTLs to exert an antigen-specific cytotoxic activity against tumour cells, whose efficacy may depend on the affinity range of TCRs and the strength of TCR-tumour antigen binding (i.e. TCR-tumour antigen binding affinity) [START_REF] Kastritis | On the binding affinity of macromolecular interactions: daring to ask why proteins interact[END_REF][START_REF] Stone | T-cell receptor binding affinities and kinetics: impact on T-cell activity and specificity[END_REF].

The presence of tumour cells expressing non-self antigens triggers the clonal expansion of CTLs with matching TCRs. Thereupon, CTL numbers are kept under control by self-regulation mechanisms [START_REF] Garrod | Dissecting T cell contraction in vivo using a genetically encoded reporter of apoptosis[END_REF][START_REF] Yu | Recent advances in CD8+ regulatory T cell research[END_REF][START_REF] Min | Spontaneous T cell proliferation: a physiologic process to create and maintain homeostatic balance and diversity of the immune system[END_REF][START_REF] Stockinger | The concept of space and competition in immune regulation[END_REF][START_REF] Troy | Cutting edge: homeostatic proliferation of peripheral T lymphocytes is regulated by clonal competition[END_REF], which play a key role in the prevention of autoimmunity. Furthermore, epigenetic and genetic processes inducing stochastic and heritable changes in the antigen expression profiles of tumour cells foster dynamical intercellular variability in the expression levels of tumour antigens [START_REF] Campoli | HLA antigen changes in malignant cells: epigenetic mechanisms and biologic significance[END_REF][START_REF] Sigalotti | Intratumor heterogeneity of cancer/testis antigens expression in human cutaneous melanoma is methylation-regulated and functionally reverted by 5-Aza-2'-deoxycytidine[END_REF][START_REF] Urosevic | Expression of melanomaassociated antigens in melanoma cell cultures[END_REF]. Due to limitations posed by self-regulation mechanisms upon the numbers of CTLs targeted against different antigens at the same tumour site, such a form of intratumour heterogeneity creates the substrate for adaptation of tumour cells to the antigen-specific cytotoxic activity of CTLs and triggers adaptive changes in the repertoire of CTLs. This results in coevolutionary dynamics whereby CTLs dynamically sculpt the antigenic distribution of tumour cells, and tumour cells concurrently reshape the repertoire of CTLs [START_REF] Phillips | Immunology taught by Darwin[END_REF].

The observation that the numbers of CD8+ and CD3+ T lymphocytes at the tumour site correlate with prognosis in different types of cancer led to the development of the 'immunoscore' as a prognostic marker in cancer patients [START_REF] Angell | From the immune contexture to the Immunoscore: the role of prognostic and predictive immune markers in cancer[END_REF][START_REF] Galon | Type, density, and location of immune cells within human colorectal tumors predict clinical outcome[END_REF][START_REF] Galon | Immunoscore and Immunoprofiling in cancer: an update from the melanoma and immunotherapy bridge 2015[END_REF][START_REF] Galon | Approaches to treat immune hot, altered and cold tumours with combination immunotherapies[END_REF]. The immunoscore provides a score that increases with the density of CD8+ and CD3+ T lymphocytes both in the centre and at the periphery of the tumour. A possible tumour classification based on the immunoscore has been proposed in [START_REF] Galon | Type, density, and location of immune cells within human colorectal tumors predict clinical outcome[END_REF], where tumours with a high immunoscore are classified as 'hot tumours', tumours with an intermediate immunoscore are classified as 'altered tumours', and tumours with a low immunoscore are classified as 'cold tumours'.

Mathematical modelling background Mathematical modelling can contribute to biomedical research on the immune response to cancer by supporting experimental results with a theoretical basis, bringing new perspectives on extant empirical data, and informing new experiments [START_REF] Altrock | The mathematics of cancer: integrating quantitative models[END_REF][START_REF] Biselli | Organs on chip approach: a tool to evaluate cancer-immune cells interactions[END_REF][START_REF] Eftimie | Interactions between the immune system and cancer: a brief review of non-spatial mathematical models[END_REF][START_REF] Eladdadi | Mathematical models of tumor-immune system dynamics[END_REF][START_REF] Handel | Simulation modelling for immunologists[END_REF][START_REF] Makaryan | Modeling immune cell behavior across scales in cancer[END_REF][START_REF] Wilkie | A review of mathematical models of cancer-immune interactions in the context of tumor dormancy[END_REF]. In particular, different aspects of the coevolutionary dynamics between immune cells and tumour cells have been studied, under the assumption of spatiallyhomogeneous mixing between cells, using a number of deterministic continuum models formulated as ordinary differential equations [START_REF] Almuallem | Oncolytic viral therapies and the delicate balance between virus-macrophage-tumour interactions: A mathematical approach[END_REF][START_REF] Aguadé-Gorgorió | Tumour neoantigen heterogeneity thresholds provide a time window for combination immunotherapy[END_REF][START_REF] Balachandran | Identification of unique neoantigen qualities in long-term survivors of pancreatic cancer[END_REF][START_REF] Cattani | Metamodeling the learning-hiding competition between tumours and the immune system: a kinematic approach[END_REF][START_REF] Onofrio | Metamodeling tumor-immune system interaction, tumor evasion and immunotherapy[END_REF][START_REF] De Pillis | A validated mathematical model of cell-mediated immune response to tumor growth[END_REF][START_REF] Frascoli | A dynamical model of tumour immunotherapy[END_REF][START_REF] Griffiths | Circulating immune cell phenotype dynamics reflect the strength of tumor-immune cell interactions in patients during immunotherapy[END_REF][START_REF] Konstorum | Addressing current challenges in cancer immunotherapy with mathematical and computational modelling[END_REF][START_REF] Kuznetsov | Nonlinear dynamics of immunogenic tumors: parameter estimation and global bifurcation analysis[END_REF][START_REF] Kuznetsov | Modeling Tumor Regrowth and Immunotherapy[END_REF][START_REF] Lin Erickson | A preliminary mathematical model of skin dendritic cell trafficking and induction of T cell immunity[END_REF][START_REF] Łuksza | A neoantigen fitness model predicts tumour response to checkpoint blockade immunotherapy[END_REF][START_REF] Mayer | Regulation of T cell expansion by antigen presentation dynamics[END_REF][START_REF] Takayanagi | A mathematical analysis of the interactions between immunogenic tumor cells and cytotoxic T lymphocytes[END_REF][START_REF] Walker | From concept to clinic: Mathematically informed immunotherapy[END_REF][START_REF] Wilkie | Mathematical models of immune-induced cancer dormancy and the emergence of immune evasion[END_REF], integro-differential equations [START_REF] Bellomo | Modelling and mathematical problems related to tumor evolution and its interaction with the immune system[END_REF][START_REF] Delitala | Recognition and learning in a mathematical model for immune response against cancer[END_REF][START_REF] Kolev | Mathematical modeling of the competition between acquired immunity and cancer[END_REF][START_REF] Lorenzi | Mathematical model reveals how regulating the three phases of T-cell response could counteract immune evasion[END_REF] and partial differential equations [START_REF] Atsou | A size and space structured model describing interactions of tumor cells with immune cells reveals cancer persistent equilibrium states in tumorigenesis[END_REF][START_REF] Al-Tameemi | Evasion of tumours from the control of the immune system: consequences of brief encounters[END_REF].

Although more amenable to analytical and numerical approaches, which allow for a more in-depth theoretical understanding of the underlying cellular dynamics, deterministic continuum models make it more difficult to incorporate the finer details of the coevolution between tumour cells and CTLs. Moreover, they cannot easily capture adaptive phenomena that are driven by stochastic effects in the evolutionary paths of single cells. Hence, one ideally wants to derive them as the appropriate limit of stochastic discrete models (i.e. individual-based models) of the coevolutionary dynamics between tumour cells and CTLs. These individual-based models track the dynamics of single cells, thus permitting the representation of single-cell-scale mechanisms, and account for possible stochastic intercellular variability in the evolutionary trajectories of individual cells. Integrating the results of computational simulations of stochastic discrete models with analytical and numerical results of their deterministic continuum counterparts makes it possible to clearly identify the validity domain of such results, thus leading to more robust biological insight.

Contents of the article

In this vein, we develop an individual-based model for the coevolutionary dynamics between tumour cells and CTLs in a well-mixed system (i.e. spatial interactions are not incorporated into the model). While being aware of the fact that a variety of different cells and molecules take part in this process, here we focus on interactions involving tumour cells and CTLs only. Every cell is viewed as an individual agent whose phenotypic state is modelled by a discrete variable. For tumour cells this variable represents a parameterisation of the antigen expression profiles, while for CTLs it represents a parameterisation of TCRs.

The model takes into account the effects of the following biological processes: proliferation and death of tumour cells and CTLs; heritable, spontaneous phenotypic changes of tumour cells resulting in variation of antigenic expression; antigen-driven expansion of CTLs (i.e. in situ clonal expansion following antigen recognition); death of tumour cells due to antigen-specific cytotoxic activity of CTLs. These processes are incorporated into the model through a set of rules that correspond to a discrete-time branching random walk on the space of phenotypic states [START_REF] Chisholm | Evolutionary dynamics of phenotype-structured populations: from individual-level mechanisms to populationlevel consequences[END_REF][START_REF] Hughes | Random walks and random environments: random walks[END_REF].

We show that a generalised version of the mathematical model of immune competition presented in [START_REF] Lorenzi | Mathematical model reveals how regulating the three phases of T-cell response could counteract immune evasion[END_REF] can be formally obtained as the deterministic continuum limit of the individualbased model presented here. This continuum model comprises a non-local partial differential equation (PDE) for the phenotype distribution of tumour cells coupled with an integrodifferential equation (IDE) for the phenotype distribution of CTLs, and shares some similarities with non-local predator-prey models such as those considered, for instance, in [START_REF] Delitala | Evolutionary branching patterns in predator-prey structured populations[END_REF][START_REF] Génieys | Adaptive dynamics: modelling Darwin's divergence principle[END_REF][START_REF] Segal | Pattern formation in a model of competing populations with nonlocal interactions[END_REF][START_REF] Tian | Nonlocal interaction driven pattern formation in a preypredator model[END_REF]. In addition to the biological phenomena incorporated into the model considered in [START_REF] Lorenzi | Mathematical model reveals how regulating the three phases of T-cell response could counteract immune evasion[END_REF], the deterministic continuum counterpart of the individual-based model developed here takes also into account the effect of changes in antigen expression profiles of tumour cells and more general forms of competitive feedback mechanisms regulating the growth of the numbers of tumour cells and CTLs.

The biologically relevant homogeneous steady-state solutions of the continuum model equations are found. Linear-stability analysis of these steady-state solutions is then carried out in order to identify possible conditions on the model parameters that may lead to different outcomes of immune competition and to the emergence of patterns of phenotypic coevolution between tumour cells and CTLs. We report on computational results of the individual-based model, and show that there is a good agreement between them and analytical and numerical results of the continuum model. Moreover, we explore possible scenarios in which differences between the outputs of the individual-based and continuum models may emerge. The results obtained disentangle the role of different cell parameters in the coevolutionary dynamics between tumour cells and CTLs.

Structure of the article

In Section 2, the individual-based model is introduced. In Section 3, the deterministic continuum counterpart of this model is presented (a formal derivation is provided in Appendix A). In Section 4, the homogeneous steady-state solutions of the continuum model equations are identified and their linear stability is investigated. In Section 5, computational results of the individual-based model are discussed and integrated with numerical solutions of the continuum model. In Section 6, key biological implications of the main findings of this study are summarised and directions for future research are outlined.

Individual-based model

We model the coevolutionary dynamics between a population of tumour cells and a population of CTLs in a well-mixed system. The population of tumour cells is labelled by the letter C, while the population of CTLs is labelled by the letter T . Building on the modelling approach developed in [START_REF] Delitala | A mathematical model for immune and autoimmune response mediated by T-cells[END_REF][START_REF] Delitala | Recognition and learning in a mathematical model for immune response against cancer[END_REF][START_REF] Lorenzi | Mathematical model reveals how regulating the three phases of T-cell response could counteract immune evasion[END_REF], at any time t ∈ [0, t f ] ⊂ R + , the phenotypic state of every tumour cell is modelled by a variable u ∈ I, where I := [-L, L] ⊂ R is the closure of the set I := (-L, L) ⊂ R with L > 0, and the phenotypic state of every CTL is modelled by a variable v ∈ I. We make the following modelling assumptions.

Assumption 2.1. The variable u represents a parameterisation of the antigen expression profiles of tumour cells, while the variable v represents a parameterisation of the target antigens of the TCRs. As a result, CTLs in the phenotypic state v will be primarily capable of recognising tumour cells in the phenotypic state u = v. Assumption 2.2. Tumour cells will have higher antigenic similarity if their phenotypic states are modelled by closer value of u. Hence, depending on the range of TCR affinity, CTLs in the phenotypic state v = u may also be capable of recognising tumour cells in phenotypic states which are sufficiently close to u. Assumption 2.3. Tumour cells in similar phenotypic states (i.e. phenotypic states that are modelled by sufficiently close values of u) will occupy similar niches and, therefore, a form of intra-population competition (i.e. clonal competition) will occur between them. Moreover, self-regulation mechanisms act on CTLs in similar phenotypic states (i.e. phenotypic states that are modelled by sufficiently close values of v). Hence, a form of intra-population competition will occur between these cells as well.

Building upon the ideas presented in [START_REF] Ardaševa | Comparative study between discrete and continuum models for the evolution of competing phenotype-structured cell populations in dynamical environments[END_REF][START_REF] Chisholm | Evolutionary dynamics of phenotype-structured populations: from individual-level mechanisms to populationlevel consequences[END_REF][START_REF] Stace | Discrete and continuum phenotype-structured models for the evolution of cancer cell populations under chemotherapy[END_REF], we model each cell as an agent that occupies a position on a lattice, which represents the space of phenotypic states. We discretise the time variable and the phenotypic states, respectively, as

t h = hτ ∈ [0, t f ], u i = iχ ∈ I and v j = jχ ∈ I, h ∈ N 0 , i, j ∈ Z,
where τ ∈ R + * and χ ∈ R + * are the time-and phenotype-step, respectively. We introduce the dependent variables N h C i ∈ N 0 and N h T j ∈ N 0 to represent, respectively, the number of tumour cells on lattice site i (i.e. in the phenotypic state u i ) and the number of CTLs on lattice site j (i.e. in the phenotypic state v j ) at time-step h. The population density functions of tumour cells and CTLs (i.e. the phenotype distributions of the two cell populations) are defined, respectively, as

n h C i ≡ n C (u i , t h ) := N h C i χ and n h T j ≡ n T (v j , t h ) := N h T j χ , (2.1)
while the total numbers of tumour cells and CTLs (i.e. the sizes of the cell populations C and T ) are defined, respectively, as

ρ h C ≡ ρ C (t h ) := i N h C i and ρ h T ≡ ρ T (t h ) := j N h T j . (2.2)
In the mathematical framework of our model, the function

I h ≡ I(t h ) := ρ T (t h ) ρ C (t h ) (2.
3) provides a possible simplified measure of the immune score at the h th time-step in the wellmixed system considered here. In particular, abstracting from the immune-score based classification of tumours recalled in Section 1, throughout the article we will refer to situations in which the average value of I, i.e. the quantity

I = τ t f h I h , (2.4)
is smaller than 1 or at least about one order of magnitude larger than 1 as 'cold tumour-like scenarios' and 'hot tumour-like scenarios', respectively, whereas the remaining situations will be classified as 'altered tumour-like scenarios'.

As mentioned earlier, we focus on a biological scenario whereby: (i) cells in the two populations divide and die due to intra-population competition (i.e. clonal competition amongst tumour cells and self-regulation of CTLs); (ii) tumour cells undergo heritable, spontaneous phenotypic changes which result in variation of antigen expression profiles; (iii) CTLs undergo antigen-driven expansion (i.e. in situ clonal expansion following antigen recognition); (iv) tumour cells die due to the antigen-specific cytotoxic activity of CTLs. These biological phenomena are incorporated into the model via the modelling strategies described in the following subsections.

Mathematical modelling of cell division and death due to intra-population competition

We assume that a dividing cell is replaced by two identical cells that inherit the phenotypic state of the parent cell (i.e. the progenies are placed on the same lattice site as their parent), while a dying cell is removed from the population.

At every time-step h, we allow tumour cells and CTLs to undergo cell division at rates α C > 0 and α T > 0, respectively. Moreover, on the basis of Assumption 2.3, at every time-step h, we allow tumour cells in the phenotypic state u i and CTLs in the phenotypic state v j to die due to intra-population competition at rates µ C K h C i and µ T K h T j , respectively, where µ C , µ T > 0 and

K h C i ≡ K C (u i , t h ) := k g(u i , u k ; θ C ) N h C k , K h T j ≡ K T (v j , t h ) := k g(v j , v k ; θ T ) N h T k . (2.5)
The function g is defined as follows g(x, y; ξ) := (2.7)

   1 |L ξ (x)| if |y -x| ≤ ξ 0 if |y -x| > ξ, for (
The quantity K h C i defined via (2.5)-(2.7) represents the number of tumour cells whose phenotypic states are modelled by values of the variable u k at a distance smaller than or equal to θ C from u i , rescaled to |L θ C (u i )|. Similarly, the quantity K h T j defined via (2.5)-(2.7) represents the number of CTLs whose phenotypic states are modelled by values of the variable v k at a distance smaller than or equal to θ T from v j , rescaled to |L θ T (v j )|. Hence, the inverse of the parameter 0 < θ C ≤ |I| (i.e. 1/θ C ) provides a measure of the level of selectivity of clonal competition amongst tumour cells and the inverse of the parameter 0 < θ T ≤ |I| (i.e. 1/θ T ) provides a measure of the level of selectivity of self-regulation mechanisms acting on CTLs (i.e. the smaller θ C and θ T , the higher the corresponding levels of selectivity). Furthermore, the parameters µ C and µ T represent the rates of death of tumour cells and CTLs due to these forms of intra-population competition.

Mathematical modelling of phenotypic changes in tumour cells

Building on the modelling strategies proposed in [START_REF] Ardaševa | Comparative study between discrete and continuum models for the evolution of competing phenotype-structured cell populations in dynamical environments[END_REF][START_REF] Chisholm | Evolutionary dynamics of phenotype-structured populations: from individual-level mechanisms to populationlevel consequences[END_REF][START_REF] Stace | Discrete and continuum phenotype-structured models for the evolution of cancer cell populations under chemotherapy[END_REF], we account for spontaneous, heritable phenotypic changes which result in variation of antigen expression profiles by allowing tumour cells to update their phenotypic states according to a random walk. More precisely, between the time-steps h and h + 1, every tumour cell enters a new phenotypic state, with probability 0 < λ C < 1, or remains in its current phenotypic state, with probability 1 -λ C . When a tumour cell in the phenotypic state u i undergoes a phenotypic change, it enters into either the phenotypic state u i-1 or the phenotypic state u i+1 with probability λ C /2. This models the fact that phenotypic changes occur randomly due to non-genetic instability, rather than being induced by selective pressures [START_REF] Huang | Genetic and non-genetic instability in tumor progression: link between the fitness landscape and the epigenetic landscape of cancer cells[END_REF]. No-flux boundary conditions are implemented by aborting any attempted phenotypic variation of a tumour cell if it requires moving into a phenotypic state outside the interval I.

Mathematical modelling of tumour-immune competition

Similarly to cell division, we assume that a CTL undergoing antigen-driven expansion is replaced by two identical cells that inherit the phenotypic state of the parent cell. Moreover, similarly to cell death due to intra-population competition, we assume that a tumour cell dying due to the antigen-specific cytotoxic activity of CTLs is removed from the population.

On the basis of Assumptions 2.1 and 2.2, at every time-step h we allow CTLs in the phenotypic state v j to undergo antigen-driven expansion at rate ζ T γ J h T j , while tumour cells in the phenotypic state u i will die due to antigen-specific cytotoxic activity of CTLs at rate

ζ C γ J h C i . Here, ζ C , ζ T , γ > 0 and J h T j ≡ J T (v j , t h ) := i g(v j , u i ; η) N h C i , J h C i ≡ J C (u i , t h ) := j g(u i , v j ; η) N h T j , (2.8)
where the function g is defined via (2.6) and (2.7). The quantity J h T j defined via (2.6)-(2.8) represents the number of tumour cells whose phenotypic states are modelled by values of the variable u i at a distance smaller than or equal to η from v j , rescaled to |L η (v j )|. Similarly, the quantity J h C i defined via (2.6)-(2.8) represents the number of CTLs in phenotypic states which are modelled by values of the variable v j at a distance smaller than or equal to η from u i , rescaled to |L η (u i )|. Hence, the parameter 0 < η ≤ |I| provides a measure of the affinity range of TCRs. In more detail, this parameter determines the range of tumour antigens each CTL can recognise: large values of η correspond to a CTL population that is able to eliminate tumour cells expressing a large spectrum of antigens, whereas low values of η correspond to scenarios where CTLs can only recognise tumour cells with a more specific antigenic expression. Furthermore, the parameter γ provides a measure of the TCR-tumour antigen binding affinity (i.e. the strength of TCR-tumour antigen binding). Previous experimental studies have highlighted the role played by TCR-tumour antigen binding affinity in anti-tumour immune response and its link with immune efficiency [START_REF] June | Car t cell immunotherapy for human cancer[END_REF][START_REF] Gerdemann | Cytotoxic t lymphocytes simultaneously targeting multiple tumor-associated antigens to treat ebv negative lymphoma[END_REF][START_REF] Tan | T cell receptor binding affinity governs the functional profile of cancer-specific cd8+ t cells[END_REF]. Finally, the parameter ζ T represents the rate at which a CTL undergoing antigen-driven expansion divides (i.e. the rate of cell division corresponding to in situ clonal expansion), and the parameter ζ C represents the rate at which a tumour cell can die due to the antigen-specific cytotoxic activity of a CTL.

Computational implementation of the individual-based model

Numerical simulations of the individual-based model are performed in MATLAB. At each timestep, we follow the procedures described hereafter to simulate phenotypic variation, cell division and death of tumour cells, and cell division and death of CTLs. The random numbers r 1 , r 2 and r 3 mentioned below are real numbers drawn from the standard uniform distribution on the interval (0, 1) using the built-in function RAND.

Computational implementation of spontaneous, heritable phenotypic changes of tumour cells

For each cell in population C, a random number, r 1 , is generated and used to determine whether the cell undergoes a phenotypic change (i.e. 0 < r 1 < λ C ) or not (i.e. λ C ≤ r 1 < 1). If the cell undergoes a phenotypic change, then a second random number, r 2 , is generated. If 0 < r 2 < 1/2, then the cell moves into the phenotypic state to the left of its current state (i.e. a cell in the phenotypic state u i will move into the phenotypic state u i-1 = u i -χ), whereas if 1/2 ≤ r 2 < 1 then the cell moves into the phenotypic state to the right of its current state (i.e. a cell in the phenotypic state u i will move into the phenotypic state u i+1 = u i +χ). As mentioned earlier, no-flux boundary conditions are implemented by aborting attempted phenotypic changes that would move a cell into a phenotypic state outside the interval I.

Computational implementation of cell division and death of tumour cells and CTLs

For each population, the number of cells in each phenotypic state is counted. The quantities K C and K T are computed via (2.5) and the quantities J C and J T are computed via (2.8), and the following definitions are used to calculate the probabilities of cell division, death and quiescence (i.e. no division nor death) for every phenotypic state of cells in populations C and T , respectively,

P b C := τ α C , P d C := τ µ C K h C i + ζ C γ J h C i , P q C := 1 -P b C + P d C (2.9)
and

P b T := τ α T + ζ T γ J h T j , P d T := τ µ T K h T j , P q T := 1 -P b T + P d T . (2.10)
Notice that we are implicitly assuming that the time-step τ is sufficiently small that 0 < P k C < 1 and 0 < P k T < 1 for all k ∈ {b, d, q}. For each cell, a random number, r 3 , is generated and the cells' fate is determined by comparing this number with the probabilities of division, death and quiescence corresponding to the phenotypic state of the cell. In more detail, for a cell in population C: if 0 < r 3 < P d C then the cell is considered dead and is removed from the population; if P d C ≤ r 3 < P d C + P b C then the cell undergoes division and an identical daughter cell is created; whereas if P d C + P b C ≤ r 3 < 1 then the cell remains quiescent (i.e. does not divide nor die). Similarly, for a cell in population T : if 0 < r 3 < P d T then the cell is considered dead and is removed from the population; if P d T ≤ r 3 < P d T + P b T then the cell undergoes division and an identical daughter cell is created; whereas if P d T + P b T ≤ r 3 < 1 then the cell remains quiescent.

Corresponding deterministic continuum model

In the case where cell dynamics are governed by the rules described in Sections 2.1-2.3, between time-steps h and h + 1 a tumour cell in the phenotypic state u i may divide, die or remain quiescent (i.e. not divide nor die) with probabilities defined via (2.9), while a CTL in the phenotypic state v j may divide, die or remain quiescent with probabilities defined via (2.10). Hence, recalling that between time-steps h and h + 1 a tumour cell in the phenotypic state u i may also enter into either of the phenotypic states u i-1 and u i+1 with probabilities λ C /2, the principle of mass balance gives the following system of coupled difference equations for the population densities n h C i and n h T j :

         n h+1 C i = 2 P b C + P q C λ C 2 n h C i+1 + n h C i-1 + (1 -λ C ) n h C i , (u i , t h ) ∈ I × (0, t f ], n h+1 T j = 2 P b T + P q T n h T i , (v j , t h ) ∈ I × (0, t f ]. (3.1)
The difference equation (3.1) 1 for n h C i is subject to no-flux boundary conditions due to the fact that, as mentioned in Section 2.2, any attempted phenotypic variation of a tumour cell is aborted if it requires moving into a phenotypic state outside the interval I.

Starting from the system of coupled difference equations (3.1), letting the time-step τ → 0 and the phenotype-step χ → 0 in such a way that

λ C χ 2 2τ → β C with 0 < β C < ∞, (3.2)
where the parameter β C is the rate of spontaneous, heritable phenotypic changes of tumour cells, using the method employed in [START_REF] Ardaševa | Comparative study between discrete and continuum models for the evolution of competing phenotype-structured cell populations in dynamical environments[END_REF][START_REF] Chisholm | Evolutionary dynamics of phenotype-structured populations: from individual-level mechanisms to populationlevel consequences[END_REF][START_REF] Stace | Discrete and continuum phenotype-structured models for the evolution of cancer cell populations under chemotherapy[END_REF], it is possible to formally show (see Appendix A) that the deterministic continuum counterpart of the stochastic discrete model comprises the following PDE-IDE system for the cell population density functions n C (u, t) and n T (v, t)

                     ∂ t n C -β C ∂ 2 uu n C = α C -µ C K C (u, t) -ζ C γ J C (u, t) n C , (u, t) ∈ I × (0, t f ], ∂ t n T = α T -µ T K T (v, t) + ζ T γ J T (v, t) n T , (v, t) ∈ I × (0, t f ], J C (u, t) := I g(u, v; η) n T (v, t) dv, K C (u, t) := I g(u, w; θ C ) n C (w, t) dw, J T (v, t) := I g(v, u; η) n C (u, t)du, K T (v, t) := I g(v, w; θ T ) n T (w, t) dw, (3.3)
with I = (-L, L). Here, the function g is defined via (2.6) and (2.7), and the non-local PDE (3.3) 1 for n C is subject to the following no-flux boundary conditions

∂ u n C (-L, t) = 0 and ∂ u n C (L, t) = 0 for all t ∈ (0, t f ]. (3.4)
We remark that linear diffusion operators like the one in the PDE (3.3) 1 , which are the deterministic, continuum counterparts of underlying random walks over the space of phenotypic states, have been widely used to model the effect of heritable, spontaneous phenotypic changes in cell populations -see, for instance, the review article [START_REF] Chisholm | Cell population heterogeneity and evolution towards drug resistance in cancer: biological and mathematical assessment, theoretical treatment optimisation[END_REF] and references therein.

Steady-state and linear-stability analyses of the continuum model equations

In this section, we first identify the biologically relevant homogeneous steady-state solutions of the continuum model equations. Then, we carry out linear-stability analysis to: (i) determine conditions that may lead to the eradication of tumour cells by CTLs or to the coexistence between the two cell populations, and (ii) identify sufficient conditions for the emergence of patterns of phenotypic coevolution between tumour cells and CTLs.

Biologically relevant steady-state solutions

A biologically relevant steady-state solution of the PDE-IDE system (3.3) subject to the boundary conditions (3.4) is given by a pair of real, non-negative functions n * C (u) and n * T (v) that satisfy the following system

   -β C ∂ 2 uu n * C = α C -µ C K * C (u) -ζ C γ J * C (u) n * C , u ∈ I, α T -µ T K * T (v) + ζ T γ J * T (v) n * T = 0, v ∈ I, (4.1)
where I = (-L, L), with (4.1) 1 subject to the boundary conditions

∂ u n * C (-L) = 0 and ∂ u n * C (L) = 0. (4.2)
In the system (4.1),

J * C (u) := I g(u, v; η) n * T (v) dv, K * C (u) := I g(u, w; θ C ) n * C (w) dw (4.3) and J * T (v) := I g(v, u; η) n * C (u) du, K * T (v) := I g(v, w; θ T ) n * T (w) dw. (4.4)
The components of homogeneous steady-state solutions satisfy the following system of equations

   α C -µ C K * C (u) -ζ C γ J * C (u) n * C = 0, u ∈ I, α T -µ T K * T (v) + ζ T γ J * T (v) n * T = 0, v ∈ I (4.5)
and are of the form

n * C (u) = ρ * C |I| ∀u ∈ I and n * T (v) = ρ * T |I| ∀v ∈ I, (4.6)
where ρ * C ≥ 0 and ρ * T ≥ 0 satisfy the following system of algebraic equations

   α C |I| -µ C ρ * C -γ C ρ * T ρ * C = 0, α T |I| -µ T ρ * T + γ T ρ * C ρ * T = 0, with γ C := ζ C γ and γ T := ζ T γ. (4.7)
The system of algebraic equations (4. In particular, since we are studying tumour-immune competition, we are interested in solutions of the system of equations (4.7) whose ρ * T component is strictly positive. There exist two solutions of this type, that is, the semitrivial solution

(ρ * C1 , ρ * T 1 ) = 0, |I|α T µ T , (4.9)
and, provided that the following condition on the model parameters is met

γ < µ T α T α C ζ C , (4.10)
the nontrivial solution

(ρ * C2 , ρ * T 2 ) = |I| (α C µ T -α T γ C ) γ T γ C + µ C µ T , |I| (α T µ C + α C γ T ) γ T γ C + µ C µ T . (4.11)
The semitrivial steady-state solution given by (4.6) and (4.9) corresponds to biological scenarios whereby tumour cells are eradicated by CTLs, while the nontrivial steady-state solution given by (4.6) and (4.11) corresponds to situations where coexistence between tumour cells and CTLs occurs. Notice that condition (4.10) indicates that lower TCR-tumour antigen binding affinity (i.e. smaller values of γ) make it more likely that tumour cells survive the cytotoxic action of CTLs, thus promoting coexistence between the two cell populations.

Linear-stability analysis

Linearising the PDE-IDE system (3.3), subject to the boundary conditions (3.4), about a steadystate of components n * C (u) and n * T (v), and using the conditions given by equations (4.1), we obtain the following PDE-IDE system for the perturbations ñC (u, t) and ñT (v, t)

     ∂ t ñC -β C ∂ 2 uu ñC = [α C -µ C K * C (u) -γ C J * C (u)] ñC -µ C KC (u, t) + γ C JC (u, t) n * C , (u, t) ∈ I × (0, t f ], ∂ t ñT = [α T -µ T K * T (v) + γ T J * T (v)] ñT -µ T KT (v, t) -γ T JT (v, t) n * T , (v, t) ∈ I × (0, t f ], (4.12)
subject to the boundary conditions

∂ u ñC (-L, t) = 0 and ∂ u ñC (L, t) = 0 for all t ∈ (0, t f ]. (4.13)
In the system (4.12), J * C (u) and K * C (u) are defined via (4.3), J * T (v) and K * T (u) are defined via (4.4), and

JC (u, t) := I g(u, v; η) ñT (v, t) dv, KC (u, t) := I g(u, w; θ C ) ñC (w, t) dw, (4.14) JT (v, t) := I g(v, u; η) ñC (u, t) du, KT (v, t) := I g(v, w; θ T ) ñT (w, t) dw. (4.15)
Due to (4.5), if the steady-state solution (n * C , n * T ) is given by (4.6) and (4.9) then the PDE-IDE system (4.12) reduces to

     ∂ t ñC -β C ∂ 2 uu ñC = α C -γ C J * C (u) ñC , (u, t) ∈ I × (0, t f ], ∂ t ñT = -µ T KT (v, t) -γ T JT (v, t) ρ * T 1 |I| , (v, t) ∈ I × (0, t f ], (4.16)
whereas if condition (4.10) is met and the steady-state solution (n * C , n * T ) is given by (4.6) and (4.11) then the PDE-IDE system (4.12) reduces to

       ∂ t ñC -β C ∂ 2 uu ñC = -µ C KC (u, t) + γ C JC (u, t) ρ * C2 |I| , (u, t) ∈ I × (0, t f ], ∂ t ñT = -µ T KT (v, t) -γ T JT (v, t) ρ * T 2 |I| , (v, t) ∈ I × (0, t f ].
(4.17)

Conditions for eradication of tumour cells by CTLs or coexistence between the two cell populations

In order to determine conditions on the model parameters that may lead to the eradication of tumour cells by CTLs or to the coexistence between the two cell populations, we study the stability of the steady-state solutions given by (4.6) and (4.9) or Substituting the ansatz (4.18) into the PDE-IDE system (4.16) and using property (4.8) along with the expression (4.9) of ρ * T 1 gives the following system of algebraic equations

     λ C = α C -γ C α T µ T C , λ T = -µ T T -γ T C α T µ T , (4.19)
which can be written in matrix form as

  α C -γ C α T µ T -λ 0 γ T α T µ T -α T -λ   C T = 0 0 .
For a non-trivial solution of system (4.19) to exist, the determinant of the above matrix must be zero. This leads to the following quadratic equation for λ

λ 2 -Bλ + C = 0 with B := α C -γ C α T µ T -α T and C := α T γ C α T µ T -α C .
Hence, the semitrivial steady-state solution given by (4.6) and (4.9) is locally asymptotically stable if the reverse of condition (4.10) holds, that is if

γ > µ T α T α C ζ C , (4.20)
since in this case B < 0 and C > 0 (i.e. Re(λ) < 0). On the other hand, performing similar calculations on the PDE-IDE system (4.17) gives the the following system of algebraic equations

       λ C = -µ C C + γ C T (α C µ T -α T γ C ) γ T γ C + µ C µ T , λ T = -µ T T -γ T C (α T µ C + α C γ T ) γ T γ C + µ C µ T , (4.21)
which can be written in matrix form as

   -µ C (α C µ T -α T γ C ) γ T γ C + µ C µ T -λ -γ C (α C µ T -α T γ C ) γ T γ C + µ C µ T γ T (α T µ C + α C γ T ) γ T γ C + µ C µ T -µ T (α T µ C + α C γ T ) γ T γ C + µ C µ T -λ    C T = 0 0 .
For a non-trivial solution of system (4.21) to exist, the determinant of the above matrix must be zero. This leads to the following quadratic equation for λ

λ 2 -Bλ + C = 0 with B := -µ C (α C µ T -α T γ C ) γ T γ C + µ C µ T + µ T (α T µ C + α C γ T ) γ T γ C + µ C µ T and C := (α C µ T -α T γ C )(α T µ C + α C γ T ) (γ T γ C + µ C µ T ) 2 [µ C µ T + γ C γ T ] .
Hence, if condition (4.10) is met, then the nontrivial steady-state solution given by (4.6) and (4.11) is locally asymptotically stable, since in this case B < 0 and C > 0 (i.e. Re(λ) < 0).

Conditions for the emergence of patterns of phenotypic coevolution between tumour cells and CTLs

In order to identify sufficient conditions for the emergence of patterns of phenotypic coevolution between tumour cells and CTLs, we study the stability of the nontrivial steady-state solution given by (4.6) and (4.11) to perturbations of the form

ñC (u, t) = C e λt ϕ k (u) and ñT (v, t) = T e λt ϕ k (v) with C , T ∈ R * , λ ∈ C. (4.22)
Here, {ϕ k } k≥1 are the eigenfunctions of the Laplace operator on I with homogeneous Neumann boundary conditions indexed by the wavenumber k, that is,

ϕ k (x) = cos (k x) with k = mπ |I| , m ∈ N, x ∈ I. (4.23)
Substituting the ansatz given by (4.22) and (4.23) into the PDE-IDE system (4.17), using the fact that

I g(x, y; ξ)ϕ k (y) dy = sin(kξ) k Ψ(x; ξ) ϕ k (x) with Ψ(x; ξ) := 2 |L ξ (x)| ,
for all x ∈ I and ξ ∈ (0, |I|], we obtain the following infinite system of algebraic equations

       λ C = -k 2 β C C -µ C sin(kθ C ) k Ψ(x; θ C ) C + γ C sin(kη) k Ψ(x; η) T ρ * C2 |I| , λ T = -µ T sin(kθ T ) k Ψ(x; θ T ) T -γ T sin(kη) k Ψ(x; η) C ρ * T 2 |I| , (4.24)
which can be written in matrix form as

    -k 2 β C -µ C sin(kθ C ) k Ψ(x; θ C ) ρ * C2 |I| -λ -γ C sin(kη) k Ψ(x; η) ρ * C2 |I| γ T sin(kη) k Ψ(x; η) ρ * T 2 |I| -µ T sin(kθ T ) k Ψ(x; θ T ) ρ * T 2 |I| -λ     C T = 0 0 .
For each x ∈ I, for a non-trivial solution of the system of algebraic equations (4.24) to exist the determinant of the above matrix must be zero. For each x ∈ I, this leads to the following quadratic equation for λ λ 2 -Bλ + C = 0 where

B ≡ B(k, x) := -k 2 β C -µ C sin(kθ C ) k ρ * C2 |I| Ψ(x; θ C ) -µ T sin(kθ T ) k ρ * T 2 |I| Ψ(x; θ T ) (4.25)
and

C ≡ C(k, x) : = k 2 β C µ T sin(kθ T ) k ρ * T 2 |I| Ψ(x; θ T ) + ρ * C2 ρ * T 2 |I| 2 γ C γ T sin(kη) k Ψ(x; η) 2 + µ C µ T sin(kθ C ) k sin(kθ T ) k Ψ(x; θ C )Ψ(x; θ T ) .
A sufficient condition for the nontrivial steady-state solution given by (4.6) and (4.11) to be driven unstable by perturbations of the form (4.22) (i.e. for patterns of phenotypic coevolution between tumour cells and CTLs to be formed) is that B > 0 and/or C < 0 so that Re(λ) > 0 for all x ∈ I. In particular, in the case where

θ C = θ T = θ, (4.26)
since k is defined via (4.23), for the condition B(k, x) > 0 to hold for all x ∈ I it suffices that

β C < 1 |I| min k∈K - sin(kθ) k ρ * C2 µ C + ρ * T 2 µ T k 2 min x∈I Ψ(x; θ),
where

K := k = mπ |I| , m ∈ N : sin(kθ) < 0 . Since, under definition (2.7), min x∈I Ψ(x; θ) = 2 max x∈I |L θ (x)| = 1 θ ,
the above condition on β C reduces to

β C < 1 |I| min k∈K - sin(kθ) kθ ρ * C2 µ C + ρ * T 2 µ T k 2 . (4.27)

Numerical simulations

In this section, we report on computational results of the individual-based model along with numerical solutions of the corresponding continuum model given by the PDE-IDE system (3.3) and subject to the boundary conditions (3.4). Simulations are integrated with the results of steady-state and linear-stability analyses of the continuum model equations presented in Section 4. In particular, we investigate the way in which the outputs of the models are affected by key parameters whose impact on the coevolutionary dynamics between tumour cells and CTLs is of particular biological interest. Such key parameters are: the TCR-tumour antigen binding affinity, γ, the level of selectivity of self-regulation mechanisms acting on CTLs, 1/θ T , the level of selectivity of clonal competition amongst tumour cells, 1/θ C , and the affinity range of TCRs, η. Moreover, we explore the existence of scenarios in which differences between the outputs produced by the two models can emerge due to effects which reduce the quality of the approximation of the individual-based model provided by the continuum model.

Set-up of numerical simulations

Without loss of generality we choose L = 1, so that I = [-1, 1] and |I| = 2, and consider a discretisation of the interval [-1, 1] consisting of 1500 points (i.e. the phenotype-step is χ ≈ 0.0013). Furthermore, we use the time-step τ = 0.05 and, unless otherwise specified, we choose the final time t f = 30 days. Building on the results of steady-state and linear-stability analyses of the continuum model equations presented in Section 4, we carry out simulations using the following initial condition for the individual-based model

n 0 C (u i ) := 10 4 (1 + a cos(A u i )), n 0 T (v j ) := 10 4 (2 + a cos(A v j )), a ≥ 0, A > 0. (5.1)
In Appendix B, we provide a detailed description of the methods employed to numerically solve the PDE-IDE system (3.3) complemented with the boundary conditions (3.4) and the continuum analogue of the initial condition (5.1), i.e. the initial condition n 0 C (u) := 10 4 (1 + a cos(A u)), n 0 T (v) := 10 4 (2 + a cos(A v)), a ≥ 0, A > 0.

(5.2)

Unless otherwise specified, we use the parameter values listed in Table 1. Here, the values of the parameters α C , α T , ζ C and ζ T are consistent with previous measurement and estimation studies on the dynamics of tumour cells and CTLs [START_REF] Christophe | A biased competition theory of cytotoxic T lymphocyte interaction with tumor nodules[END_REF][START_REF] De Pillis | Mathematical model creation for cancer chemo-immunotherapy[END_REF][START_REF] Kuznetsov | Nonlinear dynamics of immunogenic tumors: parameter estimation and global bifurcation analysis[END_REF][START_REF] Schlesinger | Coevolutionary immune system dynamics driving pathogen speciation[END_REF]. The values of the parameters µ C and µ T and the range of values of the parameters θ C and θ T are chosen so as to ensure that the equilibrium sizes and phenotype distributions of the two cell populations in isolation are biologically relevant. The range of values of the parameter η is consistent with experimental estimations of the precursor frequency of CTLs [START_REF] Blattman | Estimating the precursor frequency of naive antigen-specific CD8 T cells[END_REF], while the values of the parameter γ are consistent with those used in [START_REF] Stromberg | Robustness and fragility in immunosenescence[END_REF][START_REF] Stromberg | Diversity of T-cell responses[END_REF]. The value of the parameter λ C is taken from [START_REF] Stace | Discrete and continuum phenotype-structured models for the evolution of cancer cell populations under chemotherapy[END_REF] and corresponds to values of β C that are consistent with experimental data reported in [START_REF] Doerfler | DNA methylation: development, genetic disease and cancer[END_REF][START_REF] Duesberg | Explaining the high mutation rates of cancer cells to drug and multidrug resistance by chromosome reassortments that are catalyzed by aneuploidy[END_REF].

Main results

Eradication of tumour cells When γ is high enough so that condition (4.20) is satisfied (i.e. condition (4.10) does not hold), after initial growth, the total number of tumour cells decreases steadily over time until the tumour cell population is completely eradicated (cf. Figure 1a). This is due to the fact that, in response to a rapid growth in the size of the tumour cell population, the high TCR-tumour antigen binding affinity allows the population of CTLs to embark on rapid expansion in size that continues until CTLs have reached the critical mass required to push the population of tumour cells towards extinction. The expansion of the CTL population is followed by the transition to a contraction phase, which is characterised by a decline of the total number of CTLs to a level corresponding to the maintenance of a form of immunological memory. In fact, CTLs can persist after tumour eradication and could develop into memory T cells, thus preventing tumour outgrowth [START_REF] Prehn | Immunity to methylcholanthrene-induced sarcomas[END_REF][START_REF] Zhang | Induction of postsurgical tumor immunity and t-cell memory by a poorly immunogenic tumor[END_REF].

Hot tumour-like scenarios When γ satisfies condition (4.10) but is still sufficiently high, the total number of CTLs attains a value large enough to keep the total number of tumour cells steadily low. After initial growth, the total number of tumour cells decreases over time until it stabilises itself around a relatively small value (cf. Figure 1b). As a result, the average value of the immune score I defined via (2.3) and (2.4) is one order of magnitude larger than 1 (i.e. for the parameter values considered here I ≈ 12.7). In the framework of our model, this corresponds to the emergence of hot tumour-like scenarios.

Altered tumour-like scenarios

For intermediate values of γ that satisfy condition (4.10), after initial growth, a certain number of tumour cells and a slightly larger number of CTLs stably coexist (cf. Figure 1c). In this case, the average value of the immune score I defined via (2.3) and (2.4) is just slightly larger than 1 (i.e. for the parameter values considered here I ≈ 1.6). In the framework of our model, this corresponds to the emergence of altered tumour-like scenarios.

Cold tumour-like scenarios

For sufficiently small values of γ that satisfy condition (4.10), in the early stage of cell dynamics the total number of tumour cells overtakes the total number of CTLs, and keeps expanding until saturation (cf. Figure 1d). Accordingly, the average value of the immune score I defined via (2.3) and (2.4) is smaller than 1 (i.e. for the parameter values considered here I ≈ 0.7), which corresponds to the emergence of cold tumour-like scenarios in the framework of our model.

Robustness of numerical results

The plots in Figure 1 demonstrate that there is an excellent quantitative agreement between the results of numerical simulations of the individual-based model and numerical solutions of the corresponding continuum model. Moreover, consistently with the results of linear stability analysis of the continuum model presented in Section 4.2.1, these numerical results show that the total numbers of tumour cells and CTLs converge either to the steady-state values given by (4.9) (cf. Figure 1a), or the steady-state values given by (4.11) (cf. Figure 1b-d), depending on the fact that the choices of the model parameters are such that condition (4.20) or condition (4.10) holds, respectively. When convergence to the steady state (ρ * C2 , ρ * T 2 ) given by (4.11) occurs, in the long run, the value of the average immune score I defined via (2.4) reflects the value of the ratio ρ * T 2 /ρ * C2 . Therefore, in the framework of our tumour classification based on the average immune score I (see page 5), if condition (4.10) is met: cold tumour-like scenarios and hot tumour-like scenarios will emerge when the values of the model parameters are such that the ratio ρ * T 2 /ρ * C2 is smaller than 1 or at least about one order of magnitude larger than 1, respectively, whereas altered tumour-like scenarios will emerge in the remaining cases. This has been confirmed by the results of additional numerical simulations (results not shown). Hence, independently of the specific values of the model parameters, provided that assumption (4.10) is satisfied, cell dynamics qualitatively similar to those of Figure 1, and corresponding to hot, altered or cold tumour scenarios, will be observed depending on the value of the ratio ρ * T 2 /ρ * C2 . This testifies to the robustness of the numerical results presented here.

Patterns of phenotypic coevolution between tumour cells and CTLs: impact of the parameters θ

C and θ T Figure 2 displays the plots of the phenotype distributions of tumour cells (top panel) and CTLs (central panel) at the end of numerical simulations (i.e. close to numerical equilibrium) alongside the plots of the corresponding time evolution of the total cell numbers (bottom panel). In agreement with the analytical results presented in Section 4.2.2, when condition (4.10) is satisfied and conditions (4.26) and (4.27) are met as well, patterns of phenotypic coevolution between tumour cells and CTLs may emerge. Moreover, the top and central panels of Figure 2 show that, coherently with the shape of the function B defined via (4.25) (cf. the plots in Figure 3), smaller values of θ C and θ T correlate with the formation of more peaks in the phenotype distributions of the two cell populations. The plots in Figure 2 also demonstrate that there is a good agreement between numerical simulations of the individual-based and continuum models.

Sample temporal dynamics of such patterns are summarised by the plots in Figure 4, which show that clonal expansion leads to a rapid proliferation of CTLs that are targeted to the antigens mostly expressed by tumour cells, whereas self-regulation mechanisms induce formerly stimulated CTLs to decay. In turn, the antigen-specific cytotoxic action of CTLs causes the selection of those tumour cells that are able to escape immune recognition. As a result, immune competition induces the formation of multiple peaks in the phenotype distribution of tumour cells. This concurrently shapes the phenotype distribution of CTLs with a time shift corresponding to the time required for the CTLs to adapt to the antigenic distribution of tumour cells. The plots in Figure 4 demonstrate that there is again a good agreement between numerical simulations of the individual-based and continuum models.

Patterns of phenotypic coevolution between tumour cells and CTLs: impact of the parameter η

The results of numerical simulations summarised by the plots in Figure 5 extend the analytical results presented in Section 4.2.2 by showing that, when condition (4.10) is satisfied and η is sufficiently small, smaller values of η may induce the formation of patterns of phenotypic coevolution between tumour cells and CTLs whereby less regular multi-peaked phenotype distributions of the two cell populations emerge (cf. top and central panels of Figure 5). The temporal dynamics of such patterns are qualitatively similar to those presented in Figure 4 (results not shown). Moreover, numerical simulations indicate that smaller values of η correlate with the emergence of oscillations in the total numbers of tumour cells and CTLs, that is, CTLs undergo a succession of expansion and contraction phases that result in an alternate decay and growth of tumour cells (cf. bottom panel in Figure 5c). The plots in Figure 5 demonstrate that there is a good agreement between numerical simulations of the individual-based and continuum models.

Possible discrepancies between individual-based and continuum models

The results that have been presented so far indicate that there is a good agreement between the results of computational simulations of the individual-based model and the numerical solutions of the corresponding continuum model. However, we expect possible differences between the outputs of the two models to emerge in the presence of lower tumour cell numbers, which may lead to more pronounced demographic stochasticity, and less regular multi-peaked cell phenotype distributions, which may cause a reduction in the quality of the approximations employed in the formal derivation of the deterministic continuum model from the individual-based model. In order to investigate this, we carried out numerical simulations of the two models for choices of parameter values such that condition (4.10) holds, evolution towards relatively small tumour cell numbers occurs, and less regular cell phenotype distributions with multiple peaks emerge (see caption of Figure 6 for more details). The results obtained are summarised by the plots in Figure 6, which show that the individual-based model predicts eradication of the tumour cell population, whereas the continuum model predicts coexistence between tumour cells and 1. Initial conditions (5.1) and (5.2) with a = 1 and A = 5 were used to carry out numerical simulations. Analogous results were obtained when using different values of the parameter A (results not shown). The results from the individual-based model correspond to the average over two realisations and the related variance is displayed by the coloured areas surrounding the curves. the selectivity of antigen-specific immunity. We formally derived the deterministic continuum counterpart of such an individual-based model, and we integrated the results of numerical simulations of the two models with the results of steady-state and linear-stability analyses of the continuum model equations.

The results presented in this study shed light on the way in which different parameters shape the coevolutionary dynamics between tumour cells and CD8+ cytotoxic T lymphocytes. In particular, we demonstrated that, ceteris paribus, higher values of the TCR-tumour antigen binding affinity (i.e. the parameter γ in the model) promote the eradication of tumour cells by CTLs, while lower values facilitate the coexistence between tumour cells and CTLs. Specifically, progressively reducing the TCR-tumour antigen binding affinity brings about the emergence of: hot tumour-like scenarios, which are characterised by a large number of in situ CTLs and a low number of tumour cells, and thus represent a more fertile ground for anticancer therapeutic intervention; altered tumour-like scenarios, which reflect the intrinsic ability of the immune system to effectively mount a CTL-mediated immune response and the ability of tumour cells to partially escape such a response; cold tumour-like scenarios, which are characterised by an insufficient number of in situ CTLs and are invariably associated with poor prognosis [START_REF] Galon | Approaches to treat immune hot, altered and cold tumours with combination immunotherapies[END_REF]. This classification of tumours is also supported by experimental works showing that in situ immune reaction might be the strongest parameter influencing clinical outcome, regardless of the local tumour extension and its spread to lymph nodes [START_REF] Pagès | The essential role of the in situ immune reaction in human colorectal cancer[END_REF][START_REF] Galon | Type, density, and location of immune cells within human colorectal tumors predict clinical outcome[END_REF][START_REF] Galon | Immunoscore and Immunoprofiling in cancer: an update from the melanoma and immunotherapy bridge 2015[END_REF]. Moreover, our findings support the idea that TCR-tumour antigen binding affinity may be a good intervention target for immunotherapy that aims to turn cold or altered tumours into hot ones by enhancing CTL response. In this regard, our findings are in agreement with the conclusions of previous experimental articles indicating that a strong binding affinity of T cells to tumour antigens may play a key role in the overall immune response to the disease [START_REF] Gerdemann | Cytotoxic t lymphocytes simultaneously targeting multiple tumor-associated antigens to treat ebv negative lymphoma[END_REF].In particular, in altered tumours, increasing antigenicity, via the removal of co-inhibitory signals and/or the supply of co-stimulatory signals [START_REF] Galon | Immunoscore and Immunoprofiling in cancer: an update from the melanoma and immunotherapy bridge 2015[END_REF][START_REF] Whiteside | Emerging opportunities and challenges in cancer immunotherapy[END_REF], may enhance in situ CTLs activity, and has proven to be effective in the treatment of advanced-stage melanoma [START_REF] Wolchok | Overall survival with combined nivolumab and ipilimumab in advanced melanoma[END_REF], renal cell carcinoma [START_REF] Motzer | Nivolumab plus ipilimumab versus sunitinib in advanced renal-cell carcinoma[END_REF] and non-small cell lung cancer [START_REF] Hellmann | Nivolumab plus ipilimumab in lung cancer with a high tumor mutational burden[END_REF]. In cold tumours, a proposed approach to overcome the lack of a pre-existing immune response consists in combining a priming therapy that boosts CTL responses with the removal of co-inhibitory signals through approaches such as immune checkpoint [START_REF] Galon | Approaches to treat immune hot, altered and cold tumours with combination immunotherapies[END_REF]. The therapeutic success achieved by combining immune checkpoint therapy with chemotherapy in metastatic NSCLC has demonstrated the potential strength of this dual approach [START_REF] Gandhi | Pembrolizumab plus chemotherapy in metastatic non-small-cell lung cancer[END_REF].

Moreover, the results presented here indicate that the affinity range of TCRs (i.e. the parameter η in the model), the selectivity of clonal competition amongst tumour cells (i.e. the inverse of the parameter θ C in the model) and the selectivity of self-regulation mechanisms acting on CD8+ cytotoxic T lymphocytes (i.e. the inverse of the parameter θ T in the model) play a pivotal role in the formation of patterns of phenotypic coevolution, which create the substrate for the emergence of less regular cell phenotype distributions with multiple peaks. Such patterns are underpinned by some form of immunoediting whereby the population of CTLs evolves and continuously adapts its receptor repertoire in order to recognise and effectively eliminate tumour cells and, in turn, the antigen-specific selective pressure exerted by CTLs leads to the selection of those tumour clones that are able to evade immune recognition [START_REF] Dunn | Cancer immunoediting: from immunosurveillance to tumor escape[END_REF]. The adaptability of tumour cells and CTLs and the selective pressure they mutually exert on each other during cancer development are emerging as crucial factors in determining cancer evolutionary trajectories. This has been shown in the context of chronic lymphocytic leukemia [START_REF] Purroy | Coevolution of leukemia and host immune cells in chronic lymphocytic leukemia[END_REF] and other cancer types, as reviewed in [START_REF] George | Implications of tumor-immune coevolution on cancer evasion and optimized immunotherapy[END_REF]. Our results offer also a theoretical basis for the development of anti-cancer therapy aiming at engineering TCRs so as to shape their affinity for cancer targets [START_REF] Border | Affinityenhanced T-cell receptors for adoptive T-cell therapy targeting MAGE-A10: strategy for selection of an optimal candidate[END_REF][START_REF] Crean | Molecular Rules Underpinning Enhanced Affinity Binding of Human T Cell Receptors Engineered for Immunotherapy[END_REF][START_REF] Li | Genetically engineered T cells for cancer immunotherapy[END_REF][START_REF] Zhao | Engineered TCR-T Cell Immunotherapy in Anticancer Precision Medicine: Pros and Cons[END_REF] and adaptive therapy aiming at altering intratumour clonal competition [START_REF] Gatenby | Adaptive therapy[END_REF][START_REF] West | Towards multidrug adaptive therapy[END_REF], in order to control the coevolutionary dynamics between tumour cells and CD8+ cytotoxic T lymphocytes. In this respect, one of the best known treatment based on engineering specific TCRs is based on CAR-T cells [START_REF] Wang | Clinical manufacturing of car t cells: foundation of a promising therapy[END_REF], which confer CTLs the ability to target specific antigens. It has been demonstrated that this therapeutic strategy has several potential advantages over conventional therapies, including specificity, rapidity, high success rate and long-lasting effects [START_REF] June | Car t cell immunotherapy for human cancer[END_REF][START_REF] Gomes-Silva | Cancer immunotherapy using car-t cells: from the research bench to the assembly line[END_REF].

The good agreement between the results of numerical simulations of the individual-based and continuum models, along with the quantitative information given by (4.9) and (4.11) and the precise conditions given by (4.20) and (4.27), testifies to the robustness of the biological insight gained in this work. We also showed that possible differences between cell dynamics produced by the individual-based and continuum models can emerge under parameter settings that correspond to less regular cell phenotype distributions and more pronounced demographic stochasticity. In fact, these cause a reduction in the quality of the approximations employed in the formal derivation of the deterministic continuum model from the individualbased model (cf. Appendix A). This demonstrates the importance of integrating individualbased and continuum approaches when considering mathematical models for tumour-immune competition.

Research perspectives From a mathematical point of view, we plan to carry out a systematic investigation of the conditions on the affinity range of TCRs that may lead to the emergence of oscillations in cell numbers observed in the numerical simulations presented in this work. Moreover, from a modelling point of view, our individual-based modelling framework for the coevolutionary dynamics between tumour cells and CD8+ cytotoxic T lymphocytes, along with the formal derivation of the corresponding continuum model, can be developed further in several ways. For instance, a myriad of immunosuppressive strategies, the so-called immune checkpoints, help tumour cells acquiring features that enable them to evade immune detection, which may ultimately induce the exhaustion of CTLs in the tumour micro-environment, which impairs the immune response. The modelling approach presented here does not capture this aspect. However, exhaustion mechanisms could be incorporated into the individual-based model by, for example, allowing CTLs to enter a suppressed state (i.e. CTLs would become exhausted and thus would no longer able to eliminate tumour cells). In the continuum model, this would result in the presence of an additional loss term in the IDE (3.3) 2 along with a third equation for the dynamics of exhausted CTLs. Another track to follow to further enrich our model would be to include a spatial structure, for instance by embedding the tumour cells in the geometry of a solid tumour, and to take explicitly into account the effect of both spatial and antigen-specific interactions between tumour cells and CTLs, as similarly done in [START_REF] Kather | In silico modeling of immunotherapy and stroma-targeting therapies in human colorectal cancer[END_REF][START_REF] Macfarlane | Modelling the immune response to cancer: an individual-based approach accounting for the difference in movement between inactive and activated T cells[END_REF][START_REF] Macfarlane | A stochastic individual-based model to explore the role of spatial interactions and antigen recognition in the immune response against solid tumours[END_REF]. Including a spatial structure would make it possible, inter alia, to introduce a more precise def-inition of the immune score that incorporates the level of CTL infiltration. Furthermore, at this stage, the mathematical representation of the phenotypic state of tumour cells and CTLs employed in our modelling framework is rather abstract. This might make it difficult to carry out precise quantitative comparisons between the results of numerical simulations and experimental data. This limitation could be overcome by employing a mathematical representation of tumour antigens and TCRs similar to the one that we proposed in [START_REF] Leschiera | A mathematical model to study the impact of intra-tumour heterogeneity on anti-tumour cd8+ t cell immune response[END_REF], whereby a discrete set of tumour antigens that can be recognised by a unique repertoire of TCRs is considered. Finally, it would be interesting to incorporate explicitly into the model the effects of immunotherapeutic agents or other therapeutic agents. These are all lines of research that we will be pursuing in the future.

Appendices A Formal derivation of the continuum model

Using a method analogous to that employed in [START_REF] Ardaševa | Comparative study between discrete and continuum models for the evolution of competing phenotype-structured cell populations in dynamical environments[END_REF][START_REF] Chisholm | Evolutionary dynamics of phenotype-structured populations: from individual-level mechanisms to populationlevel consequences[END_REF][START_REF] Stace | Discrete and continuum phenotype-structured models for the evolution of cancer cell populations under chemotherapy[END_REF], we show that the PDE-IDE system (3.3) can be formally derived as the appropriate continuum limit of the individual-based model presented in this article.

Substituting definitions (2.9) of P b C and P q C into the difference equation (3.1) 1 for n h C i and definitions (2.10) of P b T and P q T into the difference equation (3.1) 2 for n h T j yields

         n h+1 C i = 1 + τ α C -τ µ C K h C i + ζ C γ J h C i λ C 2 n h C i+1 + n h C i-1 + (1 -λ C ) n h C i , n h+1 T j = 1 + τ α T + ζ T γ J h T j -τ µ T K h T j n h T i , (A.1)
where n h C i ≡ n C (u i , t h ) with (u i , t h ) ∈ I × (0, t f ] and n h T j ≡ n T (v j , t h ) with (v j , t h ) ∈ I × (0, t f ]. Using the fact that the following relations hold for τ and χ sufficiently small where the function g is defined via (2.6), the system of equations (A.1) can be formally rewritten in the approximate form If, in addition, the functions n C (u, t) and n T (v, t) are continuously differentiable with respect to the variable t, letting τ → 0 and χ → 0 in such a way that condition (3.2) is met, from the latter system of equations we formally obtain

t h ≈ t, t h+1 ≈ t + τ, u i ≈ u, u i±1 ≈ u ± χ, v j ≈ v, n h C i ≈ n C (u, t), n h+1 C i ≈ n C (u, t + τ ), n h C i±1 ≈ n C (u ± χ, t), ρ h C ≈ ρ C (t) := I n C (u, t) du,
               n C (u, t + τ ) = 1 + τ R C (K C (u, t), J C (u, t)) × × λ C 2 (n C (u + χ, t) + n C (u -χ, t)) + (1 -λ C ) n C (u, t) ,
     ∂ t n C (u, t) -β C ∂ 2 uu n C (u, t) = R C (K C , J C ) n C (u, t), (u, t) ∈ I × (0, t f ],
∂ t n T (v, t) = R T (K T , J T ) n T (v, t), (v, t) ∈ I × (0, t f ].

Substituting definitions (A.3) of R C (K C , J C ) and R T (K T , J T ) into the above system of equations gives the PDE-IDE system (3.3). Finally, the no-flux boundary conditions (3.4) follow from the fact that the attempted phenotypic variation of a tumour cell is aborted if it requires moving into a phenotypic state that does not belong to the interval I.

B Details of numerical simulations of the continuum model

To construct numerical solutions of the PDE-IDE system (3.3) subject both to the no-flux boundary conditions (3.4) and to the initial condition (5.2), we use a uniform discretisation of step ∆x = 0.0013 of the interval I = [-L, L] as the computational domain of the independent variables u and v, and a uniform discretisation of step ∆t = 0.05 of the time interval (0, t f ].

We construct numerical solutions of the non-local PDE (3.3) 1 for n C using a time-splitting approach, which is based on the idea of decomposing the original problem into simpler subproblems that are then sequentially solved at each time-step using an explicit Euler method with step ∆t. This leads to the following time-dicretisation of the PDE-IDE system (3.3) subject to the Neumann boundary conditions (3.4): where

                           n k+ 1 2 C (u) = n k C (u) + ∆t R C (K k C (u), J k C (u)) n k C (u), u ∈ [-L, L], n k+1 C (u) = n
R C (K k C , J k C ) := α C -µ C K k C -ζ C γJ k C , R T (K k T , J k T ) := α T -µ T K k T -ζ T γJ k T .
The system of equations (B.1) is numerically solved using a three-point finite difference explicit scheme for the diffusion term [START_REF] Leveque | Finite difference methods for ordinary and partial differential equations: steadystate and time-dependent problems[END_REF] and an implicit-explicit finite difference scheme for the

C Supplementary figures

  x, y; ξ) ∈ I × I × (0, |I|], (2.6) where |I| denotes the size of the interval I (i.e. |I| = 2L) and |L ξ (x)| denotes the size of the interval L ξ (x) := {y ∈ I : |y -x| ≤ ξ}, (x; ξ) ∈ I × (0, |I|].

  7) is obtained by first integrating the PDE (3.3) 1 over I and imposing the boundary conditions (3.4), then integrating the IDE (3.3) 2 over I, next substituting ansatz (4.6) into the resulting equations and equating to zero their right-hand sides, and finally using the fact that, when the function g is defined via (2.6) and (2.7), I g(x, y; ξ) dy = 1, ∀x ∈ I, ξ ∈ (0, |I|]. (4.8)

( 4 .

 4 [START_REF] Blattman | Estimating the precursor frequency of naive antigen-specific CD8 T cells[END_REF] to perturbations of the form ñC (u, t) = C e λt ∀u ∈ I and ñT (v, t) = T e λt ∀v ∈ I with C , T ∈ R * , λ ∈ C.(4.18)

Figure 1 :

 1 Figure 1: Eradication of tumour cells and emergence of hot tumour-like, altered tumourlike and cold tumour-like scenarios. Panel a. displays the plots of the time evolution of the total number of tumour cells (ρ C ) and CTLs (ρ T ) of the individual-based model (solid, coloured lines) and the continuum model (dashed, black lines) when γ is high enough that condition (4.20) is satisfied (i.e. condition (4.10) does not hold). Here, α T = 0.5/day and all the other parameters are as in Table 1 with γ = 3.5, η = 1.8 and θ C = θ T = 1.8. The grey dotted lines highlight the steady-state values of ρ C and ρ T given by (4.9). Panels b.-d. display similar plots for sufficiently large, intermediate and sufficiently small values of γ that satisfy condition (4.10) -i.e. γ = 2 (panel b.), γ = 0.3 (panel c.) and γ = 0.12 (panel d.). All the other parameters are as in Table 1 with η = 1.8 and θ C = θ T = 1.8. The grey dotted lines highlight the steady-state values of ρ C and ρ T given by (4.11). Initial conditions (5.1) and (5.2) with a = 0 were used to carry out numerical simulations. Analogous results hold when a > 0 in (5.1) and (5.2) (cf. Figure S1 in Appendix C). The results from the individual-based model correspond to the average over two realisations and the related variance is displayed by the coloured areas surrounding the curves.

Figure 2 :

 2 Figure 2: Patterns of phenotypic coevolution between tumour cells and CTLs: impact of the parameters θ C and θ T . Top panels display the plots of the population density function of tumour cells (n C ) and central panels display the plots of the population density function of CTLs (n T ) of the individual-based model (solid, coloured lines) and continuum model (dashed, black lines) at the end of numerical simulations (i.e. at t = 30) when conditions (4.10) and (4.27) are satisfied and progressively smaller values of θ C and θ T are considered -i.e. θ C = θ T = 0.5 (panels a.), θ C = θ T = 0.3 (panels b.) and θ C = θ T = 0.2 (panels c.). All the other parameters are as in Table 1 with γ = 1.5 and η = 0.7. Bottom panels display the corresponding plots of the time evolution of the total number of tumour cells (ρ C ) and CTLs (ρ T ). Initial conditions (5.1) and (5.2) with a = 1 and A = 5 were used to carry out numerical simulations. Analogous results were obtained when using different values of the parameter A (results not shown). The results from the individual-based model correspond to the average over two realisations of the underlying random walk and the related variance is displayed by the coloured areas surrounding the curves.

Figure 3 :

 3 Figure 3: Plots of B for different values of θ C and θ T . Plots of the function B(k, x) defined via (4.25) for any x ∈ argmin x∈[-1,1] Ψ(x; θ) with θ = θ C = θ T , under the parameter values used in Figure 2 -i.e. θ = θ C = θ T = 0.5 (left panel), θ = θ C = θ T = 0.3 (central panel) and θ = θ C = θ T = 0.2 (right panel).

Figure 4 :

 4 Figure 4: Sample temporal dynamics of patterns of phenotypic coevolution between tumour cells and CTLs. Top panels display the plots of the population density function of tumour cells (n C ) and bottom panels display the plots of the population density function of CTLs (n T ) of the individual-based model (solid, coloured lines) and continuum model (dashed, black lines) at five successive time instants -i.e. t = 0.4 (panels a.), t = 4 (panels b.), t = 10 (panels c.), t = 16 (panels d.), t = 30 (panels e.) -in the case where condition (4.10) is satisfied. Here, θ C = θ T = 0.3, γ = 1.5 and η = 0.7, and all the other parameters are as in Table 1. Initial conditions (5.1) and (5.2) with a = 1 and A = 5 were used to carry out numerical simulations. Analogous results were obtained when using different values of the parameter A (results not shown). The results from the individual-based model correspond to the average over two realisations and the related variance is displayed by the coloured areas surrounding the curves.

Figure 5 :

 5 Figure 5: Patterns of phenotypic coevolution between tumour cells and CTLs: impact of the parameter η. Top panels display the plots of the population density function of tumour cells (n C ) and bottom panels display the plots of the population density function of CTLs (n T ) of the individual-based model (solid, coloured lines) and continuum model (dashed, black lines) at the end of numerical simulations (i.e. at t = 30) when condition (4.10) is satisfied and progressively smaller values of η are considered -i.e. η = 1 (panels a.), η = 0.6 (panels b.) and η = 0.2 (panels c.). All the other parameters are as in Table1with γ = 1 and θ C = θ T = 0.7. Bottom panels display the corresponding plots of the time evolution of the total number of tumour cells (ρ C ) and CTLs (ρ T ). Initial conditions (5.1) and (5.2) with a = 1 and A = 5 were used to carry out numerical simulations. Analogous results were obtained when using different values of the parameter A (results not shown). The results from the individual-based model correspond to the average over five realisations of the underlying random walk and the related variance is displayed by the coloured areas surrounding the curves.

Figure 6 :

 6 Figure 6: Possible discrepancies between individual-based and continuum models. Panel a. displays the plot of the time evolution of the total number of tumour cells (ρ C ) and CTLs (ρ T ) of the individual-based model (solid, coloured lines) and the continuum model (dashed, black lines) when condition (4.10) holds, evolution towards relatively small tumour cell numbers occurs, and the parameter η is sufficiently small so that less regular multi-peaked cell phenotype distributions emerge -i.e. α T = 0.5/day, µ T = 2 × 10 -6 µl/day, γ = 1.1, η = 0.1, and all the other parameters as in Table 1 with θ C = θ T = 1.8. The plots in panels b. and c. display the corresponding population density functions of tumour cells (n C ) and CTLs (n T ) of the individual-based model (solid, coloured lines) and of the continuum model (dashed, black lines) at the end of simulations (i.e. at t = t f = 100). Initial conditions (5.1) and (5.2) with a = 1 and A = 5 were used to carry out numerical simulations. Analogous results were obtained when using different values of the parameter A (results not shown). The results from the individual-based model correspond to the average over five realisations and the related variance is displayed by the coloured areas surrounding the curves.

FundingE

  .L. has received funding from the European Research Council (ERC) under the European Union's Horizon2020 research and innovation programme (grant agreement No 740623). T.L. gratefully acknowledges support from the Italian Ministry of University and Research (MUR) through the grant "Dipartimenti di Eccellenza 2018-2022" (Project no. E11G18000350001) and the PRIN 2020 project (No. 2020JLWP23) "Integrated Mathematical Approaches to Socio-Epidemiological Dynamics" (CUP: E15F21005420006). L.A., E.L. and T.L. gratefully acknowledge support from the CNRS International Research Project "Modélisation de la biomécanique cellulaire et tissulaire" (MOCETIBI).

  J h C i ≈ J C (u, t) := I g(u, v; η) n T (v, t) dv, K h C i ≈ K C (u, t) := I g(u, w; θ C ) n C (w, t) dw, n h T j ≈ n T (v, t), n h+1 T j ≈ n T (v, t + τ ), ρ h T ≈ ρ T (t) := I n T (v, t) dv, J h T j ≈ J T (v, t) := I g(v, u; η) n C (u, t) du, K h T j ≈ K T (v, t) := I g(v, w; θ T ) n T (w, t) dw,

n

  T (v, t + τ ) = 1 + τ R T (K T (v, t), J T (v, t)) n T (v, t), (A.2) where u ∈ I, v ∈ I and t ∈ (0, t f ]. Here, R C (K C , J C ) := α C -µ C K C + ζ C γ J C , R T (K T , J T ) := α T + ζ T γ J T -µ T K T . (A.3)If the function n C (u, t) is twice continuously differentiable with respect to the variable u, for χ sufficiently small we can use the Taylor expansions n C (u ± χ, t) = n C (u,

k+ 1 2 Cuu n k+ 1 2 C

 22 (u) + ∆t β C ∂ 2 (u), u ∈ (-L, L), ∂ u n k+1 C (u) = 0, u ∈ {-L, L} n k+1 T (v) = n k T (v) + ∆t R T (K k T (v), J k T (v)) n k T (v), v ∈ [-L, L],(B.1)

Figure S1 :

 S1 Figure S1: Eradication of tumour cells and emergence of hot tumour-like, altered tumourlike and cold tumour-like scenarios: the case where a > 0. Panel a. displays the plots of the time evolution of the total number of tumour cells (ρ C ) and CTLs (ρ T ) of the individual-based model (solid, coloured lines) and the continuum model (dashed, black lines) when γ is high enough that condition (4.20) is satisfied (i.e. condition (4.10) does not hold). Here, α T = 0.5 and all the other parameters are as in Table 1 with γ = 3.5, η = 1.8 and θ C = θ T = 1.8. Panels b.-d. display similar plots for sufficiently large, intermediate and sufficiently small values of γ that satisfy condition (4.10) -i.e. γ = 2 (panel b.), γ = 0.3 (panel c.) and γ = 0.12 (panel d.). All the other parameters are as in Table 1 with η = 1.8 and θ C = θ T = 1.8. Initial conditions (5.1) and (5.2) with a = 1 and A = 5 were used to carry out numerical simulations. Analogous results were obtained when using different values of the parameter A (results not shown). The results from the individual-based model correspond to the average over two realisations and the related variance is displayed by the coloured areas surrounding the curves.

Table 1 :

 1 Parameter values used in numerical simulations and their sources

		Biological meaning	Value	Source
	αC Rate of tumour cell proliferation	1.5/day	[20]
	αT	Rate of antigen-independent CTL proliferation	5 ×10 -2 /day	[24]
	µC Rate of death of tumour cells due to clonal competition	1.5 ×10 -6 µl/day ad hoc
	µT	Rate of death of CTLs due to self-regulation mechanisms	5 ×10 -6 µl/day	ad hoc
	ζC	Killing rate of tumour cells by CTLs	5 ×10 -6 µl/day	[61]
	ζT	Rate of replication of CTLs following recognition	3 ×10 -5 µl/day	[86]
	η	Affinity range of TCRs	[0.1, 2]	[11]
	θC	Level of selectivity of clonal competition amongst tumour cells [0.1, 2]	ad hoc
	θT	Level of selectivity of self-regulation mechanisms of CTLs	[0.1, 2]	ad hoc
	γ	TCR-tumour antigen binding affinity	[0.1, 3.5]	[93, 94]
	λC Probability of phenotypic variation of tumour cells	0.01	[89]

  t) ± χ∂ u n C (u, t) + χ 2 2 ∂ 2 uu n C (u, t) + h.o.t. . (A.4) Substituting the Taylor expansions (A.4) into equation (A.2) 1 for n C (u, t+τ ), after a little algebra we find C(u, t + τ ) -n C (u, t) τ -λ C χ 2 2τ ∂ 2 uu n C (u, t) = R C (K C (u, t), J C (u, t)) n C (u, t)+ + λ C χ 2 2 R C (K C (u, t), J C (u, t)) ∂ 2 uu n C (u, t) + h.o.t. , n T (v, t + τ ) -n T (v, t) τ = R T (K T (v, t), J T (v, t)) n T (v, t).

	        	n
	       

(Tommaso Lorenzi). E.L. has received funding from the European Research Council (ERC) under the European Union's Horizon2020 research and innovation programme (grant agreement No 740623). T.L. gratefully acknowledges support of the MIUR grant "Dipartimenti di Eccellenza 2018-2022".

remaining terms [START_REF] Lorenzi | Dissecting the dynamics of epigenetic changes in phenotype-structured populations exposed to fluctuating environments[END_REF][START_REF] Lorz | Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies[END_REF], which leads to the following system of equations

Ci+1 -2n

Here,

Given the values of the parameter τ , χ and λ C of the individual-based model, the value of the parameter β C is defined so that condition (3.2) is met. The other parameter values are chosen to be coherent with those used to carry out numerical simulations of the individual-based model, which are specified in the main body of the paper.