### The Mathematical Model

Following the kinetic scheme employed by [13], in absence of immuno-editing mechanisms the interplay between tumour cells and tumour-infiltrating cytotoxic-T-lymphocytes can be modelled as shown in Figure 1 (see: [13]), where *T* denotes a tumour cell, *E* denotes an effector cell (CTL), *C* denotes the complex formed, *T*^{∗}denotes a dead tumour cell and *E*^{∗}denotes a dead effector cell. The following assumptions are made:

the complexes *C* consist of a tumour cell and a CTL forming at a rate *k*^{+} . The parameter *k*^{+} consists of the encounter rate between a tumour cell and a CTL, the probability that the CTL recognizes the tumour cell as a “non-self” entity, and also the probability that the tumour cell forms a complex with the CTLs

the break-up of complexes can lead to a situation where both the tumour cell and the CTLs are alive with a rate *k*^{−}

the break-up of complexes can lead to a situation where either the immune cell or the tumour cell survives the encounter with a rate *k*

the probability that a tumour cell is killed is *p*, and correspondingly the probability that a CTL is killed (i.e. the tumour cells survives) is (1−*p*)

Using the Law of Mass Action, this leads to the following system of differential equations describing these specific kinetic interactions:

\begin{array}{ll}{\partial}_{t}C& ={k}^{+}\text{ET}-({k}^{-}+k)C\phantom{\rule{2em}{0ex}}\\ {\partial}_{t}T& =-{k}^{+}\text{ET}+{k}^{-}C+k(1-p)C\phantom{\rule{2em}{0ex}}\\ {\partial}_{t}E& =-{k}^{+}\text{ET}+{k}^{-}C+\text{kpC}\phantom{\rule{2em}{0ex}}\end{array}

(1)

The key idea proposed in this paper which develops the work of [13], is that a proportion of the tumour cells that survive an encounter with a CTL are more resistant to any future attacks by CTLs. Consequently, the phenotypic properties of these new “enhanced” tumour cells will be different from those of the “naive” tumour cells. Specifically, we make the additional assumptions:

their probability of being killed (previously the parameter *p*) is smaller

their probability of being recognized and also of forming a complex with a CTL (embedded in the parameter *k*^{+} ) is smaller

Moreover, we shall also assume that the proliferation rate of CTLs stimulated by the presence of the complexes is also smaller. We denote the naive tumour cells by *T*_{0}(*t*) and the non-naive tumour cells by *T*_{
i
}, where *i* stands for the number of previous encounters with the CTLs. We assume that the fitness of tumour cells increases up to a maximum number of encounters *N*, implying that we consider in total 1 + *N* “classes” of tumour cells, *T*_{0},*T*_{1},…,*T*_{
N
}.

The new kinetic relationships of our model are illustrated in Figure 2and are characterized by the following new groups of parameters:

the rate of formation of complexes [*E* *T*_{
i
}]: \left(\right)close="">\n \n \n \n k\n \n \n i\n \n \n +\n \n \n \n. We assume that \left(\right)close="">\n \n \n \n k\n \n \n i\n \n \n +\n \n \n \n is constant or decreasing with index *i*, with \left(\right)close="">\n \n \n \n k\n \n \n N\n \n \n +\n \n \n \u2265\n 0\n \n;

the probability that a tumour cell of the *i*-th class is killed: *p*_{
i
}. We assume that *p*_{
i
}is decreasing with index *i*, with *p*_{
N
}≥0;

the probability of transition *T*_{
i
}→*T*_{i + 1}to the state *i*: *θ*_{
i
}. We assume that *θ*_{
i
}is increasing for 0≤*i*≤*N*−1. Since we have assumed *N* classes of tumour cells, *θ*_{
N
}=0.

As far as the temporal dynamics of the tumour cells, CTLs and complexes is concerned, once again using the Law of Mass Action, the kinetic scheme of Figure 2 can be translated into the following system of ordinary differential equations:

\begin{array}{ll}\frac{\partial {T}_{0}}{\mathrm{\partial t}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{2.56865pt}{0ex}}& =\phantom{\rule{2.56865pt}{0ex}}\phantom{\rule{0.3em}{0ex}}-{k}_{0}^{+}E{T}_{0}+{k}^{-}{C}_{0}+k(1-{\theta}_{0})(1-{p}_{0}){C}_{0}\phantom{\rule{2em}{0ex}}\\ \frac{\partial {T}_{i}}{\mathrm{\partial t}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{2.56865pt}{0ex}}& =\phantom{\rule{2.56865pt}{0ex}}\phantom{\rule{0.3em}{0ex}}-{k}_{i}^{+}E{T}_{i}+k{\theta}_{i-1}(1-{p}_{i-1}){C}_{i-1}\phantom{\rule{2em}{0ex}}\\ \phantom{\rule{2em}{0ex}}+({k}^{-}+k(1-{\theta}_{i}\left)\right(1-{p}_{i}\left)\right){C}_{i}\phantom{\rule{2em}{0ex}}\\ \frac{\partial {C}_{l}}{\mathrm{\partial t}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{2.56865pt}{0ex}}& =\phantom{\rule{2.56865pt}{0ex}}\phantom{\rule{0.3em}{0ex}}{k}_{l}^{+}E{T}_{l}-({k}^{-}+k){C}_{l}\phantom{\rule{2em}{0ex}}\\ \frac{\mathrm{\partial E}}{\mathrm{\partial t}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{2.56865pt}{0ex}}& =\phantom{\rule{2.56865pt}{0ex}}\phantom{\rule{0.3em}{0ex}}-E\left(\sum _{j=0}^{N}{k}_{j}^{+}{T}_{j}\right)+\left(\sum _{j=0}^{N}{k}^{-}{C}_{j}\right)+\left(\sum _{j=0}^{N}k{p}_{j}{C}_{j}\right)\phantom{\rule{2em}{0ex}}\end{array}

(2)

where *i*=1,…,*N*, and *l*=0,…,*N*.

However, not only the temporal but also the spatio-temporal properties of the “fitter” tumour cells are likely to be different from those of the baseline tumour cells. Namely:

the production rates of chemoattractants stimulated by a complex CTL+’non naive tumour cell’ is assumed to be smaller than that of the naive cells

since recently Vianello et al. [43] showed in an animal model that tumours produce chemicals that repels the CTLs, here we assume that those chemorepellents are produced by the non-naive cells. For the sake of the precision, Vianello’s findings showed chemoattraction for low concentrations of chemical CXCL12 emitted by tumour cells. For the sake of simplicity here we shall only consider chemorepulsion.

In the following, we provide the full equations for all variables, including the spatial components. As mentioned in the previous section, we have assumed the model of [13] as our baseline model.

#### Spatiotemporal Dynamics of the Tumour Cells

Following [13], we assume that the tumour growth may be described by a logistic law, and that the tumour cells migrate randomly. Thus, it follows that the spatio-temporal dynamics of the naive tumour cells *T*_{0} is as follows:

\begin{array}{ll}\frac{\partial {T}_{0}}{\mathrm{\partial t}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{2.77695pt}{0ex}}=& \phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{0.3em}{0ex}}\stackrel{\text{random}\phantom{\rule{1em}{0ex}}\text{motion}}{\stackrel{00A0}{{D}_{{T}_{0}}{\nabla}^{2}{T}_{0}}}+\stackrel{\text{logistic}\phantom{\rule{1em}{0ex}}\text{growth}}{\stackrel{00A0}{{r}_{1}{T}_{0}\left(1-{\beta}_{1}\sum _{j=0}^{N}{T}_{j}\right)}}\phantom{\rule{2em}{0ex}}\\ -\stackrel{\text{local}\phantom{\rule{1em}{0ex}}\text{kineticsl}}{\stackrel{00A0}{{k}_{0}^{+}E{T}_{0}+({k}^{-}+k(1-{\theta}_{0}\left)\right(1-{p}_{0}\left)\right){C}_{0}}}\phantom{\rule{2em}{0ex}}\end{array}

(3)

and the dynamics of the non-naive cells *T*_{
i
}is given by:

\begin{array}{ll}\frac{\partial {T}_{i}}{\mathrm{\partial t}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{2.77695pt}{0ex}}& =\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{0.3em}{0ex}}\stackrel{\text{random}\phantom{\rule{1em}{0ex}}\text{motion}}{\stackrel{00A0}{{D}_{{T}_{i}}{\nabla}^{2}{T}_{i}}}+\stackrel{\text{logistic}\phantom{\rule{1em}{0ex}}\text{growth}}{\stackrel{00A0}{{r}_{1}{T}_{i}\left(1-{\beta}_{1}\sum _{j=0}^{N}{T}_{j}\right)}}\phantom{\rule{2em}{0ex}}\\ \phantom{\rule{1em}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{0.3em}{0ex}}-\stackrel{\text{local}\phantom{\rule{1em}{0ex}}\text{kineticsl}}{\stackrel{00A0}{{k}_{i}^{+}E{T}_{i}+({k}^{-}+k(1-{\theta}_{i}\left)\right(1-{p}_{i}\left)\right){C}_{i}+k{\theta}_{i-1}(1-{p}_{i-1}){C}_{i-1}}}\phantom{\rule{2em}{0ex}}\end{array}

(4)

where *i*=1,…,*N*, and where *r*_{1} is the baseline exponential growth rate of the tumour (i.e. its theoretical growth rate when it is ‘small’) and *β*_{1} is the inverse of its carrying capacity (in absence of immune reactions).

#### Spatiotemporal Dynamics of the CTLs

Considering the CTLs, as in [13], both random and chemotactic motion of these cells is included. However, as previously discussed, an additional type of motility is included due to the postulated onset of “negative taxis” due to the production of a chemorepellent *ρ*by the non-naive tumour cells. This results in the following equation:

\begin{array}{ll}\phantom{\rule{-15.0pt}{0ex}}\frac{\mathrm{\partial E}}{\mathrm{\partial t}}\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{0.3em}{0ex}}=& \phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{0.3em}{0ex}}\stackrel{\text{random}\phantom{\rule{1em}{0ex}}\text{motility}}{\stackrel{00A0}{{D}_{E}{\nabla}^{2}E}}-\stackrel{\text{chemotaxis}}{\stackrel{00A0}{\chi \left(\alpha \right)\nabla .(E\nabla \alpha )}}\phantom{\rule{2em}{0ex}}\\ \phantom{\rule{0.3em}{0ex}}+\phantom{\rule{0.3em}{0ex}}\stackrel{\text{chemorepulsion}}{\stackrel{00A0}{A\left(\rho \right)\nabla .(E\nabla \rho )}}+\stackrel{\text{supply}}{\stackrel{00A0}{\text{sh}\left(x\right)}}+\stackrel{\text{proliferation}}{\stackrel{00A0}{\frac{f\sum _{j=0}^{N}{q}_{j}{C}_{j}}{g+\sum _{j=0}^{N}{T}_{j}}}}\phantom{\rule{2em}{0ex}}\\ \phantom{\rule{0.3em}{0ex}}-\phantom{\rule{0.3em}{0ex}}\stackrel{\text{decay}}{\stackrel{00A0}{{d}_{1}E}}\phantom{\rule{0.3em}{0ex}}-\phantom{\rule{0.3em}{0ex}}\stackrel{\text{local}\phantom{\rule{1em}{0ex}}\text{kineticsl}}{\stackrel{00A0}{E\left(\sum _{j=0}^{N}{k}_{j}^{+}{T}_{j}\right)\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}+\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\left(\sum _{j=0}^{N}{k}^{-}{C}_{j}\right)\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}+\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\left(\sum _{j=0}^{N}k{p}_{j}{C}_{j}\right)}}\phantom{\rule{2em}{0ex}}\end{array}

(5)

The proliferation of the CTLs, *E*, stimulated by the complexes *C*_{
j
} is embedded in the rate constant *q*_{
j
}, which, as a consequence, must be decreasing with the index *j*, with *q*_{
N
}≥0 (*f* , *g* are constant parameters). Note that in absence of immuno-editing this proliferation term reads *fC*/(*g* + *T*), and it has been has been introduced in [22, 45]. It represents the experimentally observed enhanced proliferation of CTLs in response to the tumour. This functional form is consistent with a model in which one assumes that the enhanced proliferation of CTLs is due to signals, such as released interleukins, generated by effector cells in tumour cell-CTL complexes. We note that the growth factors that are secreted by lymphocytes in complexes (e.g IL-2) act mainly in an autocrine fashion. That is to say they act on the cell from which they have been secreted and thus, in our spatial setting, their action can be adequately described by a “local” kinetic term only, without the need to incorporate any additional information concerning diffusivity. The external influx of CTLs is, for the sake of the simplicity, modelled as *sh*(*x*), where *h*(*x*) is a Heaviside function, taken to be zero over a given subregion of the domain of interest cf. [13]. In other words, we assume that there is a subdomain where lymphocytes are not naturally present and which is penetrated by CTLs only thanks to diffusion and chemotaxis.

#### Chemoattractant

The spatiotemporal dynamics of the chemoattractant *α*produced by the complexes is given by

\begin{array}{ll}\frac{\mathrm{\partial \alpha}}{\mathrm{\partial t}}\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{0.3em}{0ex}}& =\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{0.3em}{0ex}}\stackrel{\mathit{\text{diffusion}}}{\stackrel{00A0}{{D}_{2}{\nabla}^{2}\alpha}}-\stackrel{\mathit{\text{decay}}}{\stackrel{00A0}{{\delta}_{1}\alpha}}+\stackrel{\mathit{\text{production}}}{\stackrel{00A0}{\left(\sum _{j=0}^{N}{\Pi}_{j}{C}_{j}\right)}}\phantom{\rule{2em}{0ex}}\end{array}

(6)

where we assume that the production rate constant *Π*_{
i
}is decreasing with index *i*, with *Π*_{
N
}≥0, since we assume that complexes between CTLs and non-naive tumour cells are less and less able to produce such a chemoattractant.

#### Chemorepellent

Adapting the experimental findings by Vianello to our framework, we suppose that the non-naive tumour cells produce a chemical that repels the CTLs and whose concentration *ρ* is governed by the equation

\begin{array}{ll}\frac{\mathrm{\partial \rho}}{\mathrm{\partial t}}\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{0.3em}{0ex}}& =\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{0.3em}{0ex}}\stackrel{\mathit{\text{diffusion}}}{\stackrel{00A0}{{D}_{\ast}{\nabla}^{2}\rho}}-\stackrel{\mathit{\text{decay}}}{\stackrel{00A0}{{\delta}_{2}\rho}}+\stackrel{\mathit{\text{production}}}{\stackrel{00A0}{\left(\sum _{j=0}^{N}{w}_{j}{T}_{j}\right)}}\phantom{\rule{2em}{0ex}}\end{array}

(7)

where the production rate constants *w*_{
i
}are such that: *w*_{0}=0 (absence of production for naive tumour cells) and

{w}_{1}<{w}_{2}<\cdots <{w}_{N}.

(8)

#### Tumour cell-CTL Complexes

Following [13], we assume that the motility of the complexes is so small that it can be neglected:

\begin{array}{ll}\frac{\partial {C}_{l}}{\mathrm{\partial t}}\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{0.3em}{0ex}}& =\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{0.3em}{0ex}}\stackrel{\mathit{\text{local kinetic}}}{\stackrel{00A0}{{k}_{l}^{+}E{T}_{l}-({k}^{-}+k){C}_{l}}}\phantom{\rule{2em}{0ex}}\end{array}

(9)

where *l*=0,…,*N*.

#### Modeling the transition rates and probabilities

Concerning the transitions *T*_{
i
}→*T*_{i + 1}, we assume that they are a linear function of *i*:

\begin{array}{l}{\theta}_{i}={\theta}_{0}+({\theta}_{\mathit{\text{MAX}}}-{\theta}_{0})\frac{i}{N-1},\phantom{\rule{0.3em}{0ex}}i=0,\dots ,N-1\phantom{\rule{2em}{0ex}}\\ \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{1em}{0ex}}{\theta}_{N}=0,\phantom{\rule{0.3em}{0ex}}{\theta}_{\mathit{\text{MAX}}}=10{\theta}_{0}.\phantom{\rule{2em}{0ex}}\end{array}

(10)

and that their baseline value is sufficiently small: 10^{−5}≤*θ*_{0}≤10^{−3}. In other words we assume that the probability of acquiring the less immunogenic phenotype is small (in analogy with the smallness of the probability of surviving to an attack by a CTLs).

The probability *p*_{
i
} that a tumour cell of class *T*_{
i
}is lethally hit is given by:

{p}_{i}={p}_{0}+({p}_{N}-{p}_{0})\frac{i}{N},

(11)

where 0≤*p*_{
N
}<*p*_{0}. Concerning the rates \left(\right)close="">\n \n \n \n k\n \n \n i\n \n \n +\n \n \n \n, we assume either that they are constant or that they are linearly decreasing with \left(\right)close="">\n \n \n \n k\n \n \n N\n \n \n +\n \n \n =\n 0\n \n:

{k}_{i}^{+}={k}_{0}^{+}\left(1-\frac{i}{N}\right).

(12)

The production rate of the chemoattractant is also assumed to vary linearly:

{\Pi}_{i}={\Pi}_{0}\left(1-\frac{i}{N}\right)

(13)

with [16]: *Π*_{0}=20−3000*molecules* *cells*^{−1}*min*^{−1}, We suppose that the chemorepellent is produced via a mechanism of “threshold generation”, i.e. only after a sufficient number of encounters, yielding:

{w}_{i}=\left\{\begin{array}{c}0,\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\text{if}\phantom{\rule{2.77695pt}{0ex}}0\le i\le {N}_{\ast},\\ {w}_{\text{MAX}}\frac{i-{N}_{\ast}}{N-{N}_{\ast}},\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}{N}_{\ast}<i\le N\end{array}\right.

(14)

where we assumed: *w*_{
MAX
}≈*Π*_{0}(i.e. the maximum production rate of the chemorepellent is equal to the production rate of the chemoattractant in absence of immuno-editing phenomena).

### Boundary and initial conditions

We initially consider the model in a fixed 1-dimensional domain [0,*x*_{
a
}] and close the system by applying appropriate boundary and initial conditions. As far as the boundary conditions are concerned, zero-flux boundary conditions are imposed on all state variables (apart from *C*): *E*, *α*, *ρ* and *T*_{
i
}, *i*=0,…,*N*. These boundary conditions are appropriate for the tumour-immune dynamics we are considering. For example, in *BC* *L*_{1}lymphomas of the spleen tumour cells are spatially contained in the lymph tissue of the spleen, an elongated organ that, in mice, is characterized by a very strong basement membrane, which is only broken when the tumour cells switch to an “invasive phenotype”. Since here we are concerned with earlier stage dynamics of tumour cells in a dormant state evading the CTLs, it follows that the zero-flux boundary conditions are adequate for our particular model. As far as the initial conditions are concerned, we assume an initial front of naive tumour cells encountering a front of CTLs, resulting in the formation of *C*_{0}complexes. We suppose that initially there are no non-naive tumour cells and hence no complexes involving them. No chemicals are initially present in the spatial domain. These assumptions yield:

\phantom{\rule{-15.0pt}{0ex}}\begin{array}{ll}E(x,0)& =\left\{\begin{array}{ll}0,& \phantom{\rule{1em}{0ex}}\text{if}\phantom{\rule{2.56865pt}{0ex}}0\le x\le l,\\ {E}_{a}(1-\text{exp}(-1000{(x-l)}^{2}\left)\right),& \phantom{\rule{1em}{0ex}}\text{if}\phantom{\rule{2.56865pt}{0ex}}l<x\le {x}_{a}.\end{array}\right.\\ {T}_{0}(x,0)& =\left\{\begin{array}{ll}{T}_{a}(1-\text{exp}(-1000{(x-l)}^{2}\left)\right),& \phantom{\rule{1em}{0ex}}\text{if}\phantom{\rule{2.56865pt}{0ex}}0\le x\le l,\\ 0& \phantom{\rule{1em}{0ex}}\text{if}\phantom{\rule{2.56865pt}{0ex}}l<x\le {x}_{a}.\end{array}\right.\\ {C}_{0}(x,0)& =\left\{\begin{array}{ll}0,& \phantom{\rule{1em}{0ex}}\text{if}\phantom{\rule{2.56865pt}{0ex}}x\notin [l-\epsilon ,l+\epsilon ],\\ {C}_{a}\text{exp}(-1000{(x-l)}^{2}),& \phantom{\rule{1em}{0ex}}\text{if}\phantom{\rule{2.56865pt}{0ex}}x\in [l-\epsilon ,l+\epsilon ].\end{array}\right.\\ \phantom{\rule{2em}{0ex}}{T}_{i}(x,0)=0,\phantom{\rule{0.3em}{0ex}}{C}_{i}(x,0)=0,\forall x\in [0,{x}_{a}].\\ \phantom{\rule{1em}{0ex}}\phantom{\rule{2.56865pt}{0ex}}\phantom{\rule{2.56865pt}{0ex}}\phantom{\rule{2.56865pt}{0ex}}\alpha (x,0)=0,\phantom{\rule{0.3em}{0ex}}\rho (x,0)=0,\phantom{\rule{0.3em}{0ex}}\forall x\in [0,{x}_{a}].\end{array}

(15)

where

\begin{array}{l}{E}_{a}=\frac{s}{{d}_{1}},\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}{T}_{a}=\frac{1}{{\beta}_{1}},\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}{C}_{a}=\mathrm{min}({E}_{a},{T}_{a}),\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}0<\epsilon \ll 1,\phantom{\rule{2em}{0ex}}\\ i=1,\dots ,\mathrm{N.}\phantom{\rule{2em}{0ex}}\end{array}

(16)

Note that *E*_{
a
}is the baseline homogenous steady-state value of CTLs in the absence of tumour cells, and that *T*_{
a
}is the the baseline homogenous steady-state value of the naive cells in absence of CTLs, whereas *C*_{
a
} is the maximum possible density of complexes resulting from initial values of *E* and *T*.