Evasion of tumours from the control of the immune system: consequences of brief encounters

Mohannad Al-Tameemi; Mark Chaplain; Alberto d’Onofrio

doi:10.1186/1745-6150-7-31

Research
Open Access
Published: 25 September 2012

Evasion of tumours from the control of the immune system: consequences of brief encounters

Mohannad Al-Tameemi¹,
Mark Chaplain¹ &
Alberto d’Onofrio²

Biology Direct volume 7, Article number: 31 (2012) Cite this article

5294 Accesses
28 Citations
2 Altmetric
Metrics details

Abstract

Background

In this work a mathematical model describing the growth of a solid tumour in the presence of an immune system response is presented. Specifically, attention is focused on the interactions between cytotoxic T-lymphocytes (CTLs) and tumour cells in a small, avascular multicellular tumour. At this stage of the disease the CTLs and the tumour cells are considered to be in a state of dynamic equilibrium or cancer dormancy. The precise biochemical and cellular mechanisms by which CTLs can control a cancer and keep it in a dormant state are still not completely understood from a biological and immunological point of view. The mathematical model focuses on the spatio-temporal dynamics of tumour cells, immune cells, chemokines and “chemorepellents” in an immunogenic tumour. The CTLs and tumour cells are assumed to migrate and interact with each other in such a way that lymphocyte-tumour cell complexes are formed. These complexes result in either the death of the tumour cells (the normal situation) or the inactivation of the lymphocytes and consequently the survival of the tumour cells. In the latter case, we assume that each tumour cell that survives its “brief encounter” with the CTLs undergoes certain beneficial phenotypic changes.

Results

We explore the dynamics of the model under these assumptions and show that the process of immuno-evasion can arise as a consequence of these encounters. We show that the proposed mechanism not only shape the dynamics of the total number of tumor cells and of CTLs, but also the dynamics of their spatial distribution. We also briefly discuss the evolutionary features of our model, by framing them in the recent quasi-Lamarckian theories.

Conclusions

Our findings might have some interesting implication of interest for clinical practice. Indeed, immuno-editing process can be seen as an “involuntary” antagonistic process acting against immunotherapies, which aim at maintaining a tumor in a dormant state, or at suppressing it.

Reviewers

This article was reviewed by G. Bocharov (nominated by V. Kuznetsov, member of the Editorial Board of Biology Direct), M. Kimmel and A. Marciniak-Czochra.

Background

Cancer research, both experimental and theoretical, in last 12 years has been deeply influenced by the paper “The Hallmarks of Cancer” by [1]. In this paper six key aspects (“hallmarks”) of cancer development and growth were identified and examined. Their interplay with other cellular populations was mainly seen as cooperative, and thus positive for the tumour. Recently, the same authors have published an interesting follow-up paper, “Hallmarks of Cancer: The Next Generation” [2], where their description of cancer development and growth was modified and “up-dated”. In particular, among the various new topics discussed, two important new aspects were added: the role of epigenetic phenomena and the possibility of competitive interplay with the innate and adaptive immune systems. In particular, evasion from immune destruction is explicitly listed as a new “hallmark”.

Tumour cells are characterized by a large number of genetic and epigenetic events leading to the appearance of specific antigens (e.g. mutated proteins, under/over-expressed normal proteins and many others) triggering reactions by the both the innate and the adaptive immune system [3–7]. These observations have provided a theoretical basis to the empirical hypothesis of immune surveillance, i.e. that the immune system may act to eliminate tumours [8], only recently experimentally and epidemiologically confirmed [9]. Of course, the competitive interaction between tumour cells and the immune system involves a considerable number of events and molecules, and as such is extremely complex.

Moreover, to describe fully these immuno-oncological dynamics, one has to take into account a range of spatial phenomena, which are of outmost relevance in determining the dynamics of both immunogenic and non-immunogenic tumours [10–12]. In particular, the interplay between tumour cells and the immune system is strongly influenced by the spatial mobility of both tumour cells and cells of the immune system i.e. effector cells [13]. Apart from the random motion of both types of cell, a prominent role is played by chemoattraction of effector cells towards the tumour cells. Indeed, chemotactic motion of immune system cells is a hallmark of the defence of the human body against “non-self agents”, including tumours, since cells belonging to both the innate immune system (e.g. Natural Killers, Macrophages, Dendritic Cells, etc. [3]) and adaptive immune system (e.g. Cytotoxic T Lymphocites, etc. [3]) are able to reach their targets thanks to the gradients of various kinds of chemicals [3, 14], e.g. inflammation-related substances produced by tumour cells. Thus, chemotaxis is of paramount importance in the interplay between tumours and the immune system, since it influences the control of tumour growth and also the immune surveillance.

However, besides temporal and spatial non-linearities, another important point to stress is that the structure of the above-mentioned interactions is also characterized by a series of evolutionary phenomena. As is self-evident, the immune system is not able to eliminate all neoplasms. In other cases, a dynamic equilibrium may also be established, such that the tumour may survive in a so-called “dormant state” [15–18], which is undetectable (i.e. the tumour persists at a very low, undetectable level of cells but is not completely eliminated by the immune system). Until recently this was largely inferred from clinical data, but [15] have been able to show experimentally, through an ad hoc mouse model, that adaptive immunity can maintain an occult cancer in an equilibrium state. It is quite intuitive that this equilibrium can be disrupted by sudden events affecting the immune system. If disease-related impairments of the innate and adaptive immune systems, or immuno-suppressive treatments preceding organ transplantations occur, then tumour regrowth occurs [9, 19]. This has been shown both by mouse models and through epidemiological studies [9, 19].

However, there is a major class of causes of disruption of the equilibrium that is not related to immuno-suppression. Over a long period of time [9], a neoplasm may develop multiple strategies to circumvent the action of the immune system [5, 9], which may allow it to recommence growing [9, 18] into clinically apparent tumours [15], which theoretically can reach their maximum carrying capacity [18]. From an ecological point of view, we could say that the tumour has adapted to survive in a hostile environment, in which the anti-tumour immune response is activated [9, 18]. For example, the tumour may develop mechanisms to grow and spread by reducing its immunogenicity [5, 9]. In other words, the immunogenic phenotype of the tumour is influenced by the interaction with the immune system of the host. For this reason, the theory of the interactions between a tumour and the immune system has been called immuno-editing theory [9].

An impressive body of research is accumulating on immuno-evasive strategies, and a recent monograph [20] has been devoted to some aspects of this fascinating subject and to its close relationship with the effectiveness of immunotherapies. As far as the mathematical modelling of tumour and immune system interactions is concerned, there are many papers in the current literature which use deterministic models [13, 17, 18, 21–30] or stochastic models [31–35], as well as models introduced by Bellomo based on the kinetic theories of nonlinear statistical mechanics [36, 37]. The general approach of Bellomo’s theory is based on the concept of changes of activities of both the tumour cells and the effector cells of the immune system after encounters between them. As far as spatial aspects are concerned, [38, 39] developed a detailed spatio-temporal model focused on the role of macrophages. They showed that the presence of chemoattraction of macrophages towards the tumour cells implies both the onset of traveling waves and a heterogeneous spatial distribution of the tumour cells (see also [40]). Matzavinos, Chaplain and Kuznetsov proposed a spatiotemporal model of the interactions between tumour cells and cytotoxic T-lymphocytes (CTLs) [13, 16] by including the spatial motility of both tumour cells and CTLs, as well as chemotactic motion of the CTLs. They focused mainly on the role of the immune system in determining dormant states of the tumour, by showing, through a series of simulations, that a dormant state is reached, but the tumour cells are spatially distributed in an irregular pattern, which also temporally oscillates in a non-periodic fashion (see also [40]). In [18, 28], the immuno-editing phenomenon was empirically modelled by allowing the presence of slowly time-varying generic parameters in deterministic models (with time-scales significantly longer than those typical of the tumour-immune system interaction). Recently, in the framework of the above-mentioned kinetic approach, a generic model has been proposed for the learning ability of effector cells and for the hiding of tumour cells [41]. In this paper, based on the concept introduced by [9] that the immune system has the ability of “sculpting the phenotype”[9] of tumour cells (i.e. promoting the change towards less imunogenic and more resistant phenotypes), we propose a cell-centered semi-mechanistic approach aimed at describing a possible immunologically realistic kinetic mechanism through which immuno-evasion begins. Since there is strong experimental evidence that type, density and location of CTLs are predictive of the clinical outcome of some tumours, such as colorectal tumours [42], and since we are interested in the long-time dynamics, here we shall deal with the interplay of a neoplasm with CTLs. In our model, we suppose that the tumour cells that survive an attack by CTLs have a probability of acquiring (through mutations or even by epigenetic changes) a phenotype that is more resistant to future attacks by CTLs. In turn, at each new encounter with a CTL, this resistance can be increased further, and after a finite number of encounters a complete or maximal resistance to specific immunity is acquired. Moreover, specific spatial effects may be linked to the immuno-evasion of neoplasms. Indeed, recently [43] showed experimentally that tumour cells can produce chemical CXCL12 that, for large concentrations, act as chemorepellent for CTLs, whereas for low concentration it acts as a normal chemoattractor. We integrate - with some simplifications - these experimental findings in our model by permitting in the range of features defining the increasingly resistant tumour cell phenotypes an increasing ability to produce such chemorepulsive substances, whereas in a future model we shall consider the above described “nonmonotone” behavior of the taxis. These two bio-theoretical hypotheses, although new, are in line with the general schema of tumour cell escape from the immune response. Indeed, as stressed by [19], tumour cells may escape from immune control through two general mechanisms: (a) mechanisms that involve the secretion of soluble factors; (b) mechanisms that are dependent on the contact between the tumour cells and the effectors and that are aimed at reducing antigen recognition/adhesion and apoptotic resistance. Given the current experimental knowledge the above mentioned factors are primarily aimed - apart from, in many cases, their mitogenic action - at inducing the emergence of immunosuppressive networks [44]. In our present model, the factors, in line with the animal model by [43], are chemicals that repel CTLs. Finally, in the concluding remarks, we shall also discuss, from and evolutionary point of view, the differences of our model with the current evolutionary view of the immuno-editing.

Methods

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:

∂ t C = k + ET − (k − + k) C ∂ t T = − k + ET + k − C + k (1 − p) C ∂ t E = − k + ET + k − C + kpC

(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₀(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₀,T₁,…,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]: $k i +$ . We assume that $k i +$ is constant or decreasing with index i, with $k N + ≥ 0$ ;

the probability that a tumour cell of the i-th class is killed: p_i. We assume that p_iis 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 θ_iis 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:

∂ T 0 ∂t = − k 0 + E T 0 + k − C 0 + k (1 − θ 0) (1 − p 0) C 0 ∂ T i ∂t = − k i + E T i + k θ i − 1 (1 − p i − 1) C i − 1 + (k − + k (1 − θ i) (1 − p i)) C i ∂ C l ∂t = k l + E T l − (k − + k) C l ∂E ∂t = − E ∑ j = 0 N k j + T j + ∑ j = 0 N k − C j + ∑ j = 0 N k p j C j

(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₀ is as follows:

∂ T 0 ∂t = D T 0 ∇ 2 T 0 00A0 random motion + r 1 T 0 1 − β 1 ∑ j = 0 N T j 00A0 logistic growth − k 0 + E T 0 + (k − + k (1 − θ 0) (1 − p 0)) C 0 00A0 local kineticsl

(3)

and the dynamics of the non-naive cells T_iis given by:

∂ T i ∂t = D T i ∇ 2 T i 00A0 random motion + r 1 T i 1 − β 1 ∑ j = 0 N T j 00A0 logistic growth − k i + E T i + (k − + k (1 − θ i) (1 − p i)) C i + k θ i − 1 (1 − p i − 1) C i − 1 00A0 local kineticsl

(4)

where i=1,…,N, and where r₁ is the baseline exponential growth rate of the tumour (i.e. its theoretical growth rate when it is ‘small’) and β₁ 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:

∂E ∂t = D E ∇ 2 E 00A0 random motility − χ (α) ∇ . (E ∇ α) 00A0 chemotaxis + A (ρ) ∇ . (E ∇ ρ) 00A0 chemorepulsion + sh (x) 00A0 supply + f ∑ j = 0 N q j C j g + ∑ j = 0 N T j 00A0 proliferation − d 1 E 00A0 decay − E ∑ j = 0 N k j + T j + ∑ j = 0 N k − C j + ∑ j = 0 N k p j C j 00A0 local kineticsl

(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

∂α ∂t = D 2 ∇ 2 α 00A0 diffusion − δ 1 α 00A0 decay + ∑ j = 0 N Π j C j 00A0 production

(6)

where we assume that the production rate constant Π_iis 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

∂ρ ∂t = D ∗ ∇ 2 ρ 00A0 diffusion − δ 2 ρ 00A0 decay + ∑ j = 0 N w j T j 00A0 production

(7)

where the production rate constants w_iare such that: w₀=0 (absence of production for naive tumour cells) and

w 1 < w 2 < ⋯ < 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:

∂ C l ∂t = k l + E T l − (k − + k) C l 00A0 local kinetic

(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:

θ i = θ 0 + (θ MAX − θ 0) i N − 1, i = 0, …, N − 1 θ N = 0, θ MAX = 10 θ 0 .

(10)

and that their baseline value is sufficiently small: 10⁻⁵≤θ₀≤10⁻³. 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_iis lethally hit is given by:

p i = p 0 + (p N − p 0) i N,

(11)

where 0≤p_N<p₀. Concerning the rates $k i +$ , we assume either that they are constant or that they are linearly decreasing with $k N + = 0$ :

k i + = k 0 + 1 − i N .

(12)

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

Π i = Π 0 1 − i N

(13)

with [16]: Π₀=20−3000molecules cells⁻¹min⁻¹, 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 = 0, if 0 ≤ i ≤ N ∗, w MAX i − N ∗ N − N ∗, N ∗ < i ≤ N

(14)

where we assumed: w_MAX≈Π₀(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₁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₀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:

E (x, 0) = 0, if 0 ≤ x ≤ l, E a (1 − exp (− 1000 (x − l) 2)), if l < x ≤ x a . T 0 (x, 0) = T a (1 − exp (− 1000 (x − l) 2)), if 0 ≤ x ≤ l, 0 if l < x ≤ x a . C 0 (x, 0) = 0, if x ∉ [l − ε, l + ε], C a exp (− 1000 (x − l) 2), if x ∈ [l − ε, l + ε] . T i (x, 0) = 0, C i (x, 0) = 0, ∀ x ∈ [0, x a] . α (x, 0) = 0, ρ (x, 0) = 0, ∀ x ∈ [0, x a] .

(15)

where

E a = s d 1, T a = 1 β 1, C a = min (E a, T a), 0 < ε ≪ 1, i = 1, …, N.

(16)

Note that E_ais the baseline homogenous steady-state value of CTLs in the absence of tumour cells, and that T_ais 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.

Results and Discussion

We simulated our system after nondimensionalizing it as shown in the Appendix, where one can also find the numerical values of the parameters. In addition to the baseline parameter set detailed in the Appendix, in the following sections all our simulations were performed assuming the following values for key parameters associated with the encounter of CTLs and tumour cells:

θ 0 = 1 0 − 4, p N ∈ {0, 0.5, 0.75, 0.9997} .

(17)

Moreover $k 0 + = 1.3 × 1 0 − 7$ [13] and $k i +$ may be either constant or linearly decreasing, with $k N + = 0$ .

Since the average lifespan of a chimeric mouse is three years, and since we are interested in assessing the possibility (and spatio-temporal modality) of the onset of immunoevasion of a tumour, all simulations (unless stated) represent an interval of length 1100 days ≈3 years.

Spatially Homogenous Case

In this first set of simulations, we set all the spatial components of the model to zero and consider only the reaction kinetics in order to ascertain whether the primary mechanism of evasion can be purely temporal. All simulations suggest that our model, with the parameter assumptions and values we used, is able to reproduce the onset of immunoevasion in a biologically realistic time-frame. Figure 3 shows the plots of the growth of the tumour cell population over time where the killing probability at the last stage is zero, i.e. p_N=0, and $k i + = constant = 1.3 × 1 0 − 7$ . We observe that if N=4 the onset of evasion is at t≈200 days, i.e. the tumour remains dormant for 200 days, which is a long period of time for a mouse. On the contrary, if N=10 then the immunoevasion is delayed even further, with onset at t≈500 days.

Figure 4 shows the plots of the growth of the tumour cell population over time where the killing probability at the last stage is not zero but it is only halved, i.e. p_N=0.5, and as in the previous figure, $k i + = constant = 1.3 × 1 0 − 7$ . Also in this case the immunoevasion is reproduced and takes place, respectively, at t≈425 days for N=4 and at t≈950 days for N=10.

Figure 5 shows the plots of the growth of the tumour cell population over time where the killing probability at the last stage is p_N=0.75, and as in the previous figure, $k i + = constant = 1.3 × 1 0 − 7$ . These results are different from the previous two cases. Here the immunoevasion takes place in the lifespan of the mouse only for the case N=4. This suggests that in absence of changes in the parameter $k i +$ : i) the late stages T_iare the most important to determine the onset of the evasion; ii) due to the finite lifespan of chimeric mice and to the slow rate of the transitions, the immunoevasion process requires that the maximum ability of genetic or epigenetic changes in a tumour cell upon forming a complex with a CTL (embedded in the transition probability whose maximum, we recall, is at i=N−1), is reached in a small number of encounters.

Figure 6 shows the growth of the tumour cell population over time where the parameter p_N=0.75, but in this case the parameters $k i +$ are linearly decreasing with $k N + = 0$ . We notice the following differences from the previous case: i) here the onset of immunoevasion is for N=4 at t≈250 days, i.e. it is considerably accelerated; ii) there is the onset of immunoevasion (at t≈550 days) also for N=10. Thus, this simulation suggests that the role of the decrease of the probability that a tumour cell is recognized by a CTL is important for the timing of the onset of immunoevasion. Moreover, the decrease of the parameters $k i +$ alone is sufficient to induce immunoevasion, as suggested in the simulations shown in Figure 7, where p_i=constant=0.9997 and $k i +$ are linearly decreasing.

However, as shown in Figure 8, if p_N=0, then the addition of the mechanism of decreasing $k i +$ does not accelerate the onset of immunoevasion to such a degree with respect to the baseline case of constant $k i +$ shown in the previous Figure 3.

Finally, comparing the results of our simulations, in the cases where $k i +$ is decreasing we note that the maximum size of the naive tumour cells compartment is smaller than the maximum size of all the non-naive tumour cells compartments summed up:

MAX (T 0) < MAX (T 1 + ⋯ + T N) .

(18)

Moreover, Max(T₀) seems to be a decreasing function of p_N. These results might be explained as follows: the decrease of the competition between all the tumour cells and the immune system embedded in the decrease of the parameters $k i +$ , might shift the ‘internal’ competition between the naive and the non-naive tumour cells.

Spatiotemporal Model

Before we present the computational simulation results of the new model in this paper, in Figures 9 and 10 we plot the spatial distribution of tumour cells and CTLs, respectively, in the baseline case of the absence of immunoevasive mechanisms. Figure 9 shows the spatial distribution of tumour cell density within the tissue at times 100, 400, 700, and 1100 days. These results illustrate the basic spatiotemporal dynamics of the tumour cell density induced by its interplay with the distribution of CTLs i.e. a gradual transition between a front of tumour cells to a train of solitary-like travelling waves slowly invading the tissue, finally creating a spatially heterogeneous and time-changing (through irregular temporal oscillations) distribution. Similarly, Figure 10 shows the corresponding plots of the CTL density.

Figure 11 shows the spatial distribution of tumour cell density within the tissue at times 700, and 1100 days where the parameters p_N=0.75 and $k i + = constant$ . These results show that if we include the immunoevasive mechanism with a decreased value for the parameter p_N, we obtain a process that is identical to the baseline case for most of the time. Indeed, basically the left plot of this figure is identical to that of Figure9. However, after the onset of the evasion, the tumour cell density distribution changes significantly as can be seen in the left part of the domain. From these observed differences, we may say that in this modelling framework the effect of immunoevasion on the spatio-temporal dynamics is characterized by a return to a spatially homogeneous steady-state. This transition to the new spatial regimen is illustrated in more detail by the “time-slices” shown in Figure 12.

Finally, we note that the effect of chemorepulsion, which is well detectable at level of the spatial densities T_i, is no more detectable at level of the total density of all tumour cells. However, this effect is again noticeable if we consider the density of T₀ versus the density of T₁ + ⋯ + T_N(see the second plot of Figure 11).

Figures 13 and 14 show the corresponding density of CTLs in the tissue. Note that after the onset of immunoevasion, corresponding to the newly reformed invasive front of tumour cells at a high density, the density of CTLs is close to zero.

Figure 15 shows the spatial distribution of tumour cell density within the tissue at times 400, 700, and 1100 days in the case where the parameters p_N=0.75 and $k i +$ are decreasing such that $k N + = 0$ . Due to the acceleration of the immunoevasion caused by the synergy existing between the variability of p_i and $k i +$ , the plots in this figure are substantially different from those in the baseline case in Figure 9 and also with respect to the plots in Figure 11. Indeed, in this case the spatio-temporal distribution of the tumour cell density is far more regular, and by t=1100 days almost all of the tissue has been invaded by the tumour cells close to their maximum density. Moreover, here in large regions of the domain we have T₀<T₁ + ⋯ + T_N, which is the opposite of the previous case, where the naive tumour cells T₀were prevalent. Finally, Figure 15 illustrates the fact that the distributions of naive vs non-naive tumour cells are “mirror-images” of one another and they are complementary, since their sum is a homogeneous front.

In Figure 16 we show the differential effect of chemorepulsion on the various classes of tumour cells. The plots show the total number A_i(t) of cells in each classes, i.e.

A i (t) = ∫ 01 T i (x, t) dx,

(19)

over time, as well as, in the last plot, the grand-total A₁(t) + ⋯ + A_N(t). The effect of the chemorepulsion on each sub-population A_iis striking, although overall it is globally compensated (see the last plot).

Figure 17 shows the distribution of tumour cell density within the tissue at times corresponding to 700, and 1100 days respectively with parameter values p_N=0.5 and $k i + = constant = k 0 +$ . Note that in this case, although p_N=0.5, probably due to the constancy of $k i +$ , in large parts of the space the number of naive cells exceeds the rest of the classes of other tumour cells, i.e. T₀>T₁ + ⋯ + T_N. Note that at the end of the average lifespan of the mouse, the tissue is invaded to a large extent but to a lesser extent than in the case where p_N=0.5 and $k N + = 0$ . Figure 18 shows a more detailed evolution of the tumour cell density by presenting the “time-slices” from t=0 to t=1100. Figure 19 shows the corresponding distribution of CTL density.

Finally, Figure 20 shows the distribution of tumour cell density within the tissue at times corresponding to 400, 700, and 1100 days respectively when the parameters p_N=0.5 and $k N + = 0$ . This figure summarizes well the important role of the two parameters p_N and $k N +$ in shaping the spatio-temporal distribution of tumour cells. Indeed, the parameter $k N +$ appears to accelerate the onset and the velocity of propagation of the invasive front, and moreover it also differentially shapes T₀ and T₁ + ⋯ + T_N.

Conclusions

In this paper we have presented a novel mathematical model of the immune response to cancer, focusing on the specific spatio-temporal response of cytotoxic T-lymphocytes to tumour cells. We have developed and extended ideas originally formulated by [13] by proposing a possible kinetic mechanism leading to tumour evasion from the immune control. Our model is based on the key concept that a tumour cell which survives the formation of a complex with a cytotoxic T-lymphocyte can develop, with a given probability, an increased probability of surviving further attacks by CTLs. We do not specify whether this so-called ‘increased resistance’ is genetic or epigenetic. Indeed, from a kinetic point of view, this is immaterial. However, in order to experimentally validate our hypothesis, this distinction would be of paramount relevance.

Note that the model by [13] is based on a mass-action law mechanism. The reaction kinetics in a crowded cellular environment could differ from that, as studied (in absence of immunoevasion) in [17, 27, 28]. We shall investigate these more realistic settings in future works, but we think that the basic results showed here should not substantially change.

In this work we have dealt with the spatio-temporal interplay between tumours and a specific immune response from CTLs. We chose this approach because of the experimental evidence on the relevance of CTLs in determining tumour dormancy or the evasion of many important tumours such as melanomas, ovarian carcinomas and colorectal carcinomas, where the presence of infiltrating lymphocytes is a useful prognostic marker [42, 46]. However, tumour immunoevasion from dormancy is a multi-faceted phenomenon. We stress here that by no means do we think that ours is an exhaustive theoretical treatment of a such complex phenomenon.

Concerning spatial issues we also briefly mention that also in case of small dormant tumours the next generation of spatial models of tumour growth should better stress and investigate the interplay of tissue 3D geometry and the tumour vascularization with phenotypic changes in tumour cells. We have built our model based on the tumour dormancy mathematical model of [13, 16], where parameters were fitted to experimental animal (mouse) data. However, embedding the proposed evolutionary mechanism in a more complex setting, where a more detailed description of both adaptive and innate immunity is included, should lead to results qualitatively similar to those here illustrated.

Our simulations suggest that the proposed mechanism is able to mimic various dynamics of immunoevasion during the lifespan of a mouse. We have also highlighted the differential spatiotemporal contributions to evasion due, respectively, to: i) a decrease in the probability p_i of being lethally hit; ii) a decrease in the probability, embedded in $k i +$ , that a tumour cell is recognized by a CTL. In particular, our model suggests that a decrease in the parameters p_iis needed to produce evasion, which does not occur in the case where p_i remains constant at its baseline level inferred from the experimental data. However, the role of the parameters $k i +$ is important since it can greatly accelerate the simulated process. Moreover, our computational simulations also showed that the proposed mechanism can also deeply affect the spatial patterning of the tumour. In particular, our model suggests that to have a uniform invasion profile for the tumour cells necessitates also having a decrease in the recognition rate, embedded in the parameters $k i +$ . These parameters also differentially shape the spatial distribution of the various classes of tumour cells.

Concerning the possible chemorepulsion of CTLs, our computational simulation results showed that, in our biological settings, although it does not affect the spatiotemporal dynamics of the total number of tumour cells, it has a remarkable influence on the spatio-temporal distribution of the different individual classes of tumour cells. Further analysis is needed to ascertain if, with different parameters, the effect of this factor can be different, and in order to understand the behaviour in the current setting.

As far as the key ‘immuno-evasion-related’ parameters such as θ_i, p_i, and $k i +$ are concerned, we were not able to fit them with experimental data (apart, of course, from the values for p₀ and $k 0 +$ , from [13, 16]) because in the literature, to the best of our knowledge, immuno-evasion of tumours is only illustrated by means of qualitative clinical or molecular experimental findings. In particular, no immuno-evasion-related tumour growth data are available. Indeed, a complete experimental kinetic study of the adaptive evasion from tumour dormancy allowing, for example, the plotting of tumour growth curves would currently be very difficult to undertake. Thus we hope that this theoretical work may contribute to triggering such experimental investigations, which would allow us to validate our model.

From a theoretical point of view, our model, although detailed and focused on a very specific aspect of immuno-oncology, and on some very specific mechanisms, is conceptually in line with the general theories by Bellomo [36, 37, 41], who considers tumour cells and immune system effector cells as “active particles” endowed with activities and properties. Indeed, also in this paper the changes of activities of cells upon encounters between tumour cells and effector cells of the immune system are central in determining the dynamics of the system.

To the best of our knowledge, the evolutionary nature of the immuno-editing process has been studied until now under the framework of the so-called “modern synthesis”, following which the environment (in our case the immune system) is not the “causative agency” [47] but a mere selective force promoting fixation of adaptive genomic changes [47]. In the case of immunoevasion, a lowly immunogenic clone may appear spontaneously due to the large random mutation rate of tumour cells. This new phenotype is then involuntarily selected by the immune system, which kills the other phenotypes that remain strongly immunogenic. Thus the sculpting of tumour cells phenotypes [9] mentioned in the introduction, is involuntary, passive.

On the contrary, in our model the more immuno-resistant phenotypes may arise because of genetic or epigenetic causes - due to the interaction between the tumour cells and the immune system. This point of view, which is quasi-Larmackian [47], is in line with a number of recent discoveries that are leading to a new theory of “extended evolutionary synthesis” [48], which indeed investigates the impact of both genetic and epigenetic inheritance on evolutionary phenomena in order to decipher the complex interplay between genotypes, epigenotypes, phenotypes and environment. From a biological and biophysical point of view this implies also including two other key components: timescales (see, e.g., the review paper [49]) and the role of randomness. Following this extended and more modern perspective, we think that two general classes of evolutionary mechanisms might biologically underlie our kinetic model of transitions between phenotypes of tumour cells: stress-induced mutations and epigenetic switches.

In stress-induced mutagenesis [47, 49], some kind of stress induces genomic changes in a cell that, although probabilistic (e.g. in our model θ_i∈(0,1)), are triggered by the cell and its interplay with the “external world”, and that are beneficial for the cell [47, 49]. As a result, the stress-induced mutations are adaptive and beneficial, and thus they are natural candidates to explain biologically the kinetic mechanism of our mathematical model. Moreover, it is important to note that stress-induced genomic changes are a well-known and important mechanism in tumours [47, 50, 51]. Finally, Koonin hypothesized [47] that quasi-Larmarckian processes are essentially triggered by very strong signals (see e.g. Figure 3 of [47]), which, we remark, is the case for tumour-CTL interplays.

As far as the epigenetic path is concerned, it is now known that there are inheritable phenotypic changes without underlying genetic variations [48, 49], and that sometimes those changes are unusually rapid [49]. The “epigenome” is dynamic and it reflects an individual’s or a tissue’s environmental exposure during the whole of its lifespan [52]. Among the possible epigenetic changes, we mention specifically the methylation of DNA bases (“epimutations”), and also the switching between two different equilibrium states in the biochemical pathways influencing a set of inter-related phenotypes (e.g. “low immunogenicity” and “large immunogenicity”), due for example, to strong stress signals. Indeed, due to hysteresis phenomena inducing memory, also when an external stress signal is removed the system may not return to its original state. This memory may remain stable for many generations thus providing an example of epigenetic inheritance [49].

In this work we were interested essentially in the basic facts of the immune response to tumours. However, a number of immunotherapies have been proposed and also theoretically investigated (see [17, 26, 30] and references therein). We believe that both the experimental results concerning immunoevasion of tumours and the theoretical findings we have proposed here might have some implication of interest for clinical practice. Indeed, the immunoevasion mechanisms may be seen as a process counter-acting immunotherapies able to maintain a tumour in a dormant state. In general we share the opinion of [46], who stressed that recent progress in immuno-oncology has not influenced the way anti-cancer therapies are conceived and applied clinically. Indeed, we think that the immuno-editing process can be seen as an “involuntary” antagonistic process acting against immunotherapies. As far as this point is concerned, in future we will also investigate the possibility that the presence of an immunotherapy is somewhat sensed by a tumour, with a consequential acceleration of immuno-editing.

The model that was proposed here has to be understood as a detailed model at the level of the kinetics of the cellular populations of a possible mechanism that might enable tumour cells to evade from the control of adaptive immunity. The various specific (and tumour-dependent) strategies deployed by those cells in order to reach their aim, are phenomenologically described by means of the model of the dependence of the various parameters on the classes of tumour cells, as well as in the macroscopic modelling of the chemorepulsion. This is a first step in a research effort for a more complete description of tumour cell immunoevasion, which will include the detailed modelling of the biological mechanisms underlying those and other specific evasion strategies. Thus, given the complex network of interplaying between inter-cellular and intra-cellular signalling, and given the various temporal scales (from the rapid dynamics of the intracellular pathways involved, to the relatively slow growth of a tumour, up to the very slow onset of immunoevasion) as well as the spatial ones (from individual cells to visible tumours), a more detailed model will have to be multiscale. This will involve a wide array of computational tools, from those typical of computational biology and bioinformatics, to more classical analytical and numerical methods of statistical mechanics and mathematical physics.

Appendix

Non-dimensionalization

Before undertaking any computational simulations, we non-dimensionalise our model (as well as the boundary and initial conditions) by adopting the following scaling: i) space is scaled by adopting a reference value x_a equal to the size of the domain in consideration and we assume x_a=1cm; time is rescaled relative to the diffusion rate of CTLs by setting $t a = x a 2 D E$ :

x ̄ = x x a, t ̄ = t t a;

(20)

ii) cellular densities are rescaled relative to the maxima of the respective initial conditions:

T ̄ i = T i T a, i = 0, …, N Ē = E E a, C ̄ l = C l C a, l = 0, …, N,

(21)

where, we recall, E_a=C_a; iii) The concentrations of the chemoattractant and of the chemorepellent are rescaled relative to the baseline values α_a and ρ_arespectively:

α ̄ = α α a, ρ ̄ = ρ ρ a

(22)

By omitting the bars, for the sake of simplifying the notation, the proposed model becomes:

∂E ∂t = ∇ 2 E − γ (α) ∇ (E ∇ α) + ξ (ρ) ∇ (E ∇ ρ) + η 1 ∑ j = 0 N q j C j a + ∑ j = 0 N T j − E ∑ j = 0 N ψ j T j − σE + k 3 ∑ j = 0 N C j + k 4 ∑ j = 0 N p j C j + σh (x), ∂α ∂t = D α ∇ 2 α − δ α α + μ α ∑ j = 0 N Π j C j, ∂ρ ∂t = D ρ ∇ 2 ρ − δ ρ ρ + μ ρ ∑ j = 0 N w j T j, ∂ T 0 ∂t = ∇ 2 T 0 + r T 0 1 − ∑ j = 0 N T j − ϕ 0 E T 0 + k 1 C 0 + (1 − θ 0) (1 − p 0) k 2 C 0, ∂ T i ∂t = ∇ 2 T i + r T i 1 − ∑ j = 0 N T j − ϕ i E T i + k 1 C i + θ i − 1 (1 − p i − 1) k 2 C i − 1 + k 2 (1 − θ i) (1 − p i) C i, ∂ C l ∂t = ψ l E T l − λ C l,

(23)

where:

r = r 1 t a ϕ i = k i + E a t a k 1 = k − C a t a T a a = g T a λ = t a (k − + k) ψ i = k i + T a t a k 2 = k C a t a T a γ = α a t a χ η 1 = f t a T a σ = μ o t a k 3 = k − t a μ α = C a t a α a D α = D 2 t a δ α = δ 1 t a k 4 = k t a μ ρ = T a t a ρ a ξ = α a t a A D ρ = D ∗ t a δ ρ = δ 2 t a

(24)

After the non-dimensionalization, the boundary conditions become (in 1-D):

∂E ∂x (0, t) = 0, ∂E ∂x (1, t) = 0, ∂α ∂x (0, t) = 0, ∂α ∂x (1, t) = 0, ∂ρ ∂x (0, t) = 0, ∂ρ ∂x (1, t) = 0, ∂ T i ∂x (0, t) = 0, ∂ T i ∂x (1, t) = 0,

(25)

which then imply, assuming some smoothness of the solution and the form of C_iequations, that

∂ C i ∂x (0, t) = 0, ∂ C i ∂x (1, t) = 0

(26)

and the initial conditions take the form (again in 1-D):

E (x, 0) = 0, if 0 ≤ x ≤ l (1 − exp (− 1000 (x − l) 2)), if l < x ≤ 1 α (x, 0) = 0, ∀ x ∈ [0, 1] ρ (x, 0) = 0, ∀ x ∈ [0, 1] T 0 (x, 0) = (1 − exp (− 1000 (x − l) 2)), if 0 ≤ x ≤ l 0, if l < x ≤ 1 T i (x, 0) = 0, ∀ x ∈ [0, 1] C 0 (x, 0) = 0, if x ∉ [l − ε, l + ε] exp (− 1000 (x − l) 2), if x ∈ [l − ε, l + ε] C i (x, 0) = 0, ∀ x ∈ [0, 1]

(27)

where

l = 0.2, ε = 0.01 .

(28)

The chosen non-null initial conditions for the naive tumour and immune cells represent a front of tumour cells encountering a front of CTLs, resulting in the formation of CTL-tumour cell complexes. In the absence of a tumour, the homogeneous steady-state density of the CTLs is s/d₁ and therefore this is the value we have taken for the initial density E_aof CTLs in the initial conditions. Similarly, in the absence of an immune response, the homogeneous steady-state density of the tumour cells is 1/β₁ and this is what we take as the initial density of tumour cells T_ain the initial conditions. Thus, when the fronts of the two cell populations meet, the maximum density of CTL-tumour cell complexes will be min(E_a,T_a) and hence our choice for C_a. The function exp(−1000(x−l)²) was chosen to mimick the onset of sharp but continuous front.

Numerical values of the parameters

To carry out the computational simulations of the proposed model, we used the baseline parameter set reported in [13], since these were all estimated from experimental data on murine B-cell Lymphoma BC L₁, which is an animal model for the study of tumour dormancy [53]. In addition to this baseline set, we used the migration parameters proposed in [16]. Thus, the complete set of parameters, and of their meaning, is the following (TC stands for tumour Cell):

s = 1.36 × 104 da y − 1 cells c m − 1 CTL supply rate f = 0.2988 × 108 da y − 1 cells c m − 1 Complex-induced proliferatiomn rate constant g = 2.02 × 107 cells c m − 1 Constant appearing in proliferation function k 0 + = 1.3 × 1 0 − 7 da y − 1 cell s − 1 cm ‘Formation rate’ of complexes k − = 24 da y − 1 Complexes ’unbinding rate’ k = 7.2 da y − 1 Rate of lethally hit TC or CTL rate p 0 = 0.9997 Death probability of naive TC r 1 = 0.18 da y − 1 Baseline exponentil growht rate of TCs β 1 = 2 × 1 0 − 9 cell s − 1 cm Inverse of tumour carrying capacity d 1 = 0.0412 da y − 1 Baseline death rate of CTLs D E = 1 0 − 6 c m 2 da y − 1 Diffusion coefficient of CTLs χ = 1.728 × 106 c m 2 da y − 1 M − 1 Chemotaxis coefficient of CTLs D T i = 1 0 − 6 c m 2 da y − 1 Diffusion coefficient of TCs D 2 = 8 × 1 0 − 3 c m 2 da y − 1 Diffusion coefficient of chemorattractant δ 1 = 1.1551 0 − 2 day s − 1 degradation rate of the chemoattractant

(29)

Hence, from the experimental data above, the non-dimensional values of parameters becomes:

γ = 1.728 × 102 η 1 = 5.976 × 104 ψ 0 = 6.5 × 107 σ = 4.12 × 104 k 3 = 2.4 × 107 k 4 = 7.2 × 106 D α = 8 × 103 δ α = 1.115 × 104 μ α Π 0 = 104 r = 1.8 × 105 ϕ 0 = 4.29 × 104 k 1 = 1.584 × 104 k 2 = 4.752 × 103 λ = 3.12 × 107 a = 4.04 × 1 0 − 2

(30)

As far as the spatiotemporal dynamics of the chemorepellent is concerned, we assume that its diffusion coefficient and decay rate is the same as the chemokine αi.e. D_ρ=D_α, δ_ρ=δ_α, μ_ρ=μ_α, and ξ=γ.

Note that in the first table above we did not provide the value of Π₀, wherease in the second table we directly provided the adimensional value μ_αΠ₀=10⁴, in line with [13, 16].

Reviewers’ comments

The comments of the referees were reported in italics. All the three referees included minor comments on misprints or undefined parameters or other minor suggestions, which were all implemented. Thus, we only reported the comments of interest to the general readership.

Note for the referees: following your suggestions the revised version now contains an Appendix.

Reviewer #1: Prof. G. Bocharov (nominated by Dr. V. Kuznetsov, member of the Editorial Board of Biology Direct) (Institute of Numerical Mathematics, Russian Academy of Sciences, Moscow, Russian Federation)

The paper presents a theoretical study of the tumour immuno-evasion from dormancy. To this end a mathematical model of reaction-diffusion-chemotaxis type formulated with PDE is proposed. The model considers the spatio-temporal population dynamics of interactions between tumour cells and cytotoxic T cells. The key assumption of the model reflecting recent biological insights into the pathogenesis of the solid tumour growth states that the tumour cell population is heterogeneous with respect to the parameters characterizing the outcome of the interaction with cytotoxic T lymphocytes. The interaction with CTL appears to acts as a selection force shaping the microevolution of the tumour. The heterogeneity assumption enters the model via parameterized dependence of the rate of transition of tumour cells from naive to more mature states, the efficacy of CTL mediated killing and the production of chemicals acting as attractants and repellors for CTLs. Numerical simulations with the model for various parameters combinations show that the immuno-evasion can result from the phenotypic heterogeneity of the tumour cells dynamically adapting to and shaped by the anti-tumour immune response.

General remarks:

1. The issue of spatio-temporal modelling of tumour growth is an area of active research in cell population dynamics (e.g., the studies of Bertuzzi and Gandolfi). The related work on spatio-temporal modelling of tumour growth needs to be refereed to in the Introduction section.

We fully agree, and we added the required references.

Page 6: The authors make use of the Mass Action Law to describe the interaction between CTL and target cells via bilinear terms. Meantime, the reaction kinetics in crowded cell environments can differ from those assumed by the classical chemical kinetics. This issue need to be commented.

This point is very important, and we commented it in the conclusions.

Page 10: The parameter “1000” appearing in the initial functions has to be justified. Please, explain the choice of C _a .

We added a full justification of the chosen initial conditions.

Page 13: For the readership, the parameters of the model should be presented in a Table with columns specifying notation, biological meaning, units and uncertainty ranges.

We added a table with average values of the parameters, as obtained in Ref. 13 by means of least squares algorithm, and their meaning.

Page 13, line 20: Can the parameters equality assumption for the repellor and attractant be biologically justified?

Unfortunately there is a lack of experimental results on this important data, so this is a pure assumption.

Page 14, line -5: Please, elaborate more on the specific choice of the parameter values.

We added a reference for this parameter. More in general all parameters in absence of immuno-editing were taken as in references [13, 16].

Page 15-16: Figures3-8could be presented in a more compact way (e.g., as array of graphs). The same applies to Figures9-19.

We tried your suggestion, but the result was not less clear. Take also into the account that Biology Direct is purely online, thus there are no pages restrictions.

Page 18: Section 2-Dimensional Domain is based upon simulations with the chemotaxis and chemorepulsion not included in the model. The value of the results of the 2D model based simulations is not straightforward. Please, either expand the section by considering the same set of processes as in 1D case or remove it.

We removed it.

Conclusions section: It would be interesting to have the authors opinion on whether the tissue 3D geometry and its vascularisation represent important structural constraints affecting the genomic or phenotypic changes in tumour cells that need to be considered in the next generation of the spatially resolved models of tumour growth?

We fully agree, and we mentioned this important issue in the discussions.

Finally, all suggestions in your minor remarks were implemented.

Reviewer #2: Prof. M. Kimmel (Rice University, Houston, USA)

The paper concerns an extension of Chaplain group’s model (Matzavinos et al. 2004) of response of tumour cells to cytotoxic T-lymphocytes. The extension involves classes from naive through non-naive tumour cells, which have different ability to repel T-lymphocytes and other differing characteristics. The models consists of a complicated system of partial differential equations and is investigated purely by simulation.

In this first investigation we employed numerical analysis for two main reasons: first, we feel the immuno-editing phenomena as mainly a transitory effect, and as such requiring mainly numerical analysis; second, we wanted to write a manuscript aimed to theoretical biologists with biological background, such as the general readership of ’Biology Direct’. However, we also outlined a more compact version of the proposed model where the phenotypes are described by means of a continuum. For this version we soon start a full mathematical analysis.

In my opinion, stratification of tumour cells with respect to their immuno-naivete is a natural concept, which however, leads to many complications and ambiguities in model building. The paper only partially reviews the consequences of varying the hypotheses.

In the above mentioned continuum version we should be able of better exploring these important issues.

In particular, the 2-d model by the end of the paper seems superfluous, as it is not thoroughly investigated and seems to be plagued by boundary artifacts. I suggest removing it for the time being.

Following you suggestions, and the suggestions of another referee, we removed it.

I am not sure if the tumour immuno-editing can be considered a Lamarckian evolutionary process. Even if it is, I am not sure if this is important for the paper.

We respectfully disagree with you on this point. The possibility that immuno-editing might be, at least in some cases, a Lamarckian process, where the interplay between the different kind of cells is the ’driver’, is of some interest, in our opinion.

Page 2. I do not know if “extremely complex” implies “strongly non-linear” as the authors seem to claim.

We removed this observation.

Page 2.(bottom). It would be perhaps helpful to specify which cell types belong to the adaptive immune system and which to the innate one.

Done.

Page 5. “Since we are modeling a situation where immuno-evasion of the tumour cells is not considered …” Discussion of what this assumption really implies will be helpful.

We removed the above quite ambiguous phrase. Indeed, we wanted to means that the baseline model, with the considered values of its parameters, referred to a case where the tumour is not able (without immuno-editing) to evade the immune control.

Pages 9-13. Part of this material can probably go into an online Appendix.

Done.

Reviewer #3: Prof. A. Marciniak-Czochra (Heidelberg University)

The paper is devoted to modelling of spatio-temporal dynamics of interactions of cytotoxic T-lymphocytes with cells of a homogeneous avascular tumour. It takes into consideration a chemoattraction of the immune cells by the tumour and a chemorepulsion to chemicals produced by non-naive tumour cells. The model is based on the one presented in the paper “Mathematical Modelling of the spatio-temporal Response of Cytotoxic T-Lymphocyte to a Solid Tumour” (Matzavinos A, Chaplain M, Kuznetsov, Mathematical Medicine and Biology, 2004, 21:1-34) The novelty of the current work is in the assumption that the tumour cells, surviving encounter with the immune cells, undergo a functional change, which benefits their immune resistance. This resistance increases with each encounter until the maximal resistance is reached. The resulting system of coupled partial and ordinary di erential equations is solved numerically. The results are discussed in the biological context. I recommend the manuscript for publication in Biology Direct after some minor revisions related to the organisation of the paper. The issues, which should be taken into account before publication are (Note of authors:this referee listed a long series of misprints, which were all fixed. We only report the comments of interest for the reader):

P.8, the second Equation: f and g are not introduced. Although model equations are nearly identical to those in reference (Matzavinos et al, 2004), authors should motivate the term “proliferation and recruitment”.

In the revision, we fully explained the meaning of the term like fC/(g + T).

P.10: It would be helpful to have the initial conditions plotted.

Since the shape of the initial conditions are elementry, and since another referee asked to add some written comemnts on this point, we did not include the requested plots.

P.13: It should be mentioned how parameters are estimated.

The parameters were estimanted in Kuznetsov et al 1991 basically by means of mininmal least squares estimate algorithm.

Please motivate the choice of θ₀; p_N; $k 0 +$ ; how sensitive are the results to the choice of the parameters?

We have choosen a small value for θ₀because we think that the probability of acquiring the fitter phenotype is small, in analogy with the small probability of surviving to an attack by a CTL. The value of $k 0 +$ is taken from (Matzavinos et al 2004).

P. 28, Figure 12: Why (700,1100) and not (800, 1100)?

800 days would have been a too long time with respect to 700 (100 days in mice correspond approximatively to 2300 days in men).

Axis labels of the following figures are not readable: 9, 10, 11, 13, 15, 17, 20, 21.

Concerning the readability of the axis labels, it probably depend on the fact that our eps figures were converted by the publisher in jpeg figures. In our original pdf figures they seems readable. We shall carefully manage this point in case the manuscript should be accepted.

The following were comments on the revised manuscript (other comments on missprints are omitted, and were implemented):

The rescaling of equations have mistakes. For example on p. 30, λ should be $t a (k − + k)$ ,

Indeed, it was a misprint. We apologize.

while $ψ i = k i + E a T a t a / C i$ .

We had forgot to stress that C_a=E_a, so what we wrote is correct.

I find the following phrase hardly readable “Finally, we note that the predicted effect on the tumour cells spatial distribution of the induction of chemorepulsion of CTLs is not detectable if we consider the total density of all tumour cells.” Please refomulate it.

We reformulated it.

Author’s contributions

AD and MC defined the mathematical model; MAT performed the simulations; AD, MC and MAT discussed the results; AD,MC and MAT wrote the manuscript. All authors read and approved the final manuscript.

References

1.
Hanahan D, Weinberg RA: The Hallmarks of Cancer. Cell. 2000, 100: 57-70. 10.1016/S0092-8674(00)81683-9.
PubMed CAS Article Google Scholar
2.
Hanahan D, Weinberg RA: Hallmarks of Cencer: The Next Generation. Cell. 2011, 144: 646-674. 10.1016/j.cell.2011.02.013.
PubMed CAS Article Google Scholar
3.
Delves P, Martin S, Burton D, Roitt I: Roitt’s Essent Immunol. 2011, Chichester: Wiley-Blackwell
Google Scholar
4.
Kindt T, Osborne B, Goldsby R: Kuby Immunol. 2006, New York: W. H. Freeman
Google Scholar
5.
Pardoll D: Does the Immune System See Tumors as Foreign or Self?. Annu Rev Immununol. 2003, 21: 807-839. 10.1146/annurev.immunol.21.120601.141135.
CAS Article Google Scholar
6.
Swann J, Smyth M: Immune Surveillance of Tumors. J Clininical Investigations. 2007, 117: 1137-1146. 10.1172/JCI31405.
CAS Article Google Scholar
7.
Yefenov E: Innate and Adaptive Immunity in the Tumor Microenvironment. 2008, New York: Springer
Google Scholar
8.
Ehrlich P: Ueber den jetzigen stand der karzinomforschung. Nederlands Tijdschrift voor Geneeskunde. 2009, 53: 273-290.
Google Scholar
9.
Dunn G, Old L, Schreiber R: The Three Es of Cancer Immunoediting. Annu Rev Immununology. 2004, 22: 329-360. 10.1146/annurev.immunol.22.012703.104803.
CAS Article Google Scholar
10.
Adams V, Bellomo N: A Survey of Models for Tumor-Immune Syst Dyn. 1996, Boston: Birkhauser
Google Scholar
11.
Bertuzzi A, d’ Onofrio A, Fasano A, Gandolfi A: Modelling cell populations with spatial structure: steady state and treatment-induced evolution of tumour cords. Discrete and Continuous Dynamical Syst-Ser B. 2004, 4: 115-134.
Google Scholar
12.
Andasari V, Gerisch A, Lolas G, South A, Chaplain M: Mathematical modeling of cancer cell invasion of tissue: biological insight from mathematical analysis and computational simulation. J Math Biol. 2011, 63: 141-171. 10.1007/s00285-010-0369-1.
PubMed Article Google Scholar
13.
Matzavinos A, Chaplain M, Kuznetsov V: Mathematical Modelling of the spatio-temporal Response of Cytotoxic T-Lymphocyte to a Solid Tumour. Math Med Biol: A J of the IMA. 2004, 21: 1-34. 10.1093/imammb/21.1.1.
Article Google Scholar
14.
Keener J, Sneyd J: Math Physiol. 2003, New York: Springer-Verlag
Google Scholar
15.
Koebel CM, Vermi W, Swann JB, Zerafa N, Rodig SJ, Old LJ, Smyth MJ, Schreiber RD: Adaptive immunity maintains occult cancer in an equilibrium state. Nat. 2007, 450: 903-907. 10.1038/nature06309.
CAS Article Google Scholar
16.
Chaplain M, Matzavinos A: Mathematical Modelling of Spatio-temporal Phenomena In Tumour Immunology. Lecture Notes in Math. 2006, 1872: 131-183.
Article Google Scholar
17.
d’Onofrio A: A genaral Framework for Modeling Tumor-Immune System Competition and Immunotherapy: Analysis and Medical Inferences. Physica D. 2005, 208: 220-235. 10.1016/j.physd.2005.06.032.
Article Google Scholar
18.
d’Onofrio A: Tumor Evasion from Immune Control: Strategies of a MISS to Become a MASS. Chaos, Solitons and Fractal. 2007, 31: 261-268. 10.1016/j.chaos.2005.10.006.
Article Google Scholar
19.
Stewart TJ, Abrams SI: How Tumour Escape Mass Destruction. Oncogene. 2008, 27: 5894-5903. 10.1038/onc.2008.268.
PubMed CAS Article Google Scholar
20.
Gabrilovich DI, Hurwitz AA: Tumor-Induced Immune Suppression. 2008, New York: Springer
Google Scholar
21.
Stepanova N: Course of the Immune Reaction During the Development of a Malignant Tumor. Biophy. 1980, 24: 917-923.
Google Scholar
22.
Kuznetsov VA, Makalkin IA, Taylor MA, Perelson AS: Nonlinear Dynmamics of Immunogenic Tumors: Parameter Estimation and Global Bifercation Analysis. Bull Math Biol. 1994, 56 (2): 295-321.
PubMed CAS Article Google Scholar
23.
Kuznetsov V, Knott G: Modeling Tumor Regrowth and Immunotherapy. Math Comput Modell. 2001, 33: 1275-1287. 10.1016/S0895-7177(00)00314-9.
Article Google Scholar
24.
de Pillis L, Radunskaya A, Wiseman C: A Validated Mathematical Model of Cell-Mediated Immune Response to Tumor Growth. Cancer Res. 2005, 65: 7950-7958.
PubMed CAS Google Scholar
25.
Kogan Y, Forys U, Shukron O, Kronik N, Agur Z: Cellular Immunotherapy for High Grade Gliomas: Mathematical Analysis Deriving Efficacious Infusion Rates Based on Patient Requirement. SIAM J Appl Math. 2010, 70: 1953-1976. 10.1137/08073740X.
CAS Article Google Scholar
26.
Arciero JC, Jackson TL, Kirschner DE: A Mathematical Model of Tumour-Immune Evasion and siRNA Treatment. Discrete and Continuous Dynamical Syst Ser B. 2004, 4: 39-58.
Google Scholar
27.
d’Onofrio A: The Role of the Proliferation Rate of Effectors in the Tumor-Immune System Competition. Math Mododels and Methods in Appl Sci. 2006, 16: 1375-1401. 10.1142/S0218202506001571.
Article Google Scholar
28.
d’Onofrio A: Metamodeling Tumor-Immune System Interaction, Tumor Evasion and Immunotherapy. Math Comput Modell. 2008, 47: 614-637. 10.1016/j.mcm.2007.02.032.
Article Google Scholar
29.
Kronick N, Kogan Y, Elishmereni M, Halevi-Tobias K, Vuk Pavlovic S, Agur Z: Predicting Effect of Prostate Cancer Immunotherapy By Personalized Mathematical Models. PLoS One. 2010, 5: e15482-10.1371/journal.pone.0015482.
Article Google Scholar
30.
Kim PS, Lee PP, Levy D: Dynamics and Potential Impact of the Immune Response to Chronic Myelogenous Leukemia. PLoS Computuational Biol. 2008, 4: e1000095-10.1371/journal.pcbi.1000095.
Article Google Scholar
31.
Horsthemke W, Lefever R: Phase Transitions Induced By External Noise. Phys Lett A. 1977, 64: 19-21. 10.1016/0375-9601(77)90512-6.
Article Google Scholar
32.
Du LC, Mei DC: The Critical Phenomenon and the Re-entrance Phenomenon in the Anti-Tumor Model Induced By the Time Delay. Phys Lett A. 2010, 374: 3275-3279. 10.1016/j.physleta.2010.06.003.
CAS Article Google Scholar
33.
Caravagna G, d’Onofrio A, Milazzo P, Barbuti R: Antitumour Immune Surveillance Through Stochastic Oscillations. J Theor Biol. 2010, 265 (3): 336-345. 10.1016/j.jtbi.2010.05.013.
PubMed CAS Article Google Scholar
34.
d’Onofrio A: Tumour Evasion from Imune System Control as Bounded-Noise Induced Transition. Phys Rev E. 2010, 81: art. n. 021923-
Google Scholar
35.
Caravagna G, d’Onofrio A, Barbuti R: Fine-tuning anti-tumor immunotherapies via stochastic simulations. BMC Bioinf. 2012, 13 (Suppl 4): art.n. S8-
Article Google Scholar
36.
Bellomo N, Bellouquid A, Delitala M: MAthematical Topics on the Modelling Complex Multicellular Systems and Tumor Immune Cell Competition. Math Models and Methods in Appl Sci. 2004, 14: 1683-1733. 10.1142/S0218202504003799.
Article Google Scholar
37.
Bellomo N, Delitala M: From the Mathematical Kinetic, and Stochastic Game Theory for Active Particles to Modelling Mutations, Onset, Progression and Immune Competition of Cancer Cells. Phys Life Rev. 2008, 5: 183-206. 10.1016/j.plrev.2008.07.001.
Article Google Scholar
38.
Owen M, Sherratt J: Pattern Formation and Spatiotemporal Irregularity in a Model for Macrophage-Tumour Interactions. J Theor Biol. 1997, 189: 63-80. 10.1006/jtbi.1997.0494.
PubMed CAS Article Google Scholar
39.
Owen M, Sherratt J: Mathematical Modelling of Macrophage Dynamics in Tumours. Math Methods and Models in Appl Sci. 1999, 9: 513-539. 10.1142/S0218202599000270.
Article Google Scholar
40.
d’Onofrio A: Spatiotemporal effects of a possible chemorepulsion of tumor cells by immune system effectors. J Theoreoretical Biol. 2012, 296: 41-48.
Article Google Scholar
41.
Bellomo N: Modelling the Hiding–learning Dynamics in Large Living Systems. Appl Math Lett. 2010, 23: 907-911. 10.1016/j.aml.2010.04.008.
Article Google Scholar
42.
Galon J, Costes A, Sanchez-Cabo F, Kirilovsky A, Mlecnik B, Lagorce-Pagès C, Tosolini M, Camus M, Berger A, Wind P, Zinzindohoué F, Bruneval P, Cugnenc PH, Trajanoski Z, Fridman WH, Pagès F: Type, density, and location of immune cells within human colorectal tumors predict clinical outcome. Sci. 2006, 313: 1960-1964. 10.1126/science.1129139.
CAS Article Google Scholar
43.
Vianello F, Papeta N, Chen T, Kraft P, White N, Hart WK, Kircher MF, Swart E, Rhee S, Irimia D, Toner M, Weissleder R, Poznansky MC, Palù G: Murine B16 melanomas expressing high levels of the chemokine stromal-derived factor-1/CXCL12 induce tumor-specific T cell chemorepulsion and escape from immune control. J Immunol. 2006, 176: 2902-2914.
PubMed CAS Article Google Scholar
44.
Kim R, Emi M, Tanabe K: Cancer Immunnoediting from Immune Surveillance to Immune Escape. Immunol. 2007, 121: 1-14. 10.1111/j.1365-2567.2007.02587.x.
CAS Article Google Scholar
45.
Kuznetsov V: A mathematical model for the interaction between cytotoxic lymphocytes and tumor cells. Analysis of the growth, stabilization and regression of the B cell lymphoma in mice chimeric with respect to the major histocompatibility complex. Biomed Sci. 1991, 2: 465-476.
PubMed CAS Google Scholar
46.
Zitvogel L, Tesniere A, Kroemer G: Cancer Despite Immunosurveillance: Immunoselection and Immunosubversion. Nat Rev Immunol. 2006, 6: 715-727. 10.1038/nri1936.
PubMed CAS Article Google Scholar
47.
Koonin EV, Wolf YI: Is evolution Darwinian and/or Lamarckian?. Biol Direct. 2009, 4: 42-10.1186/1745-6150-4-42.
PubMed PubMed Central Article Google Scholar
48.
Danchin E, Charmantier A, Champagne F, Mesoudi A, Pujol B, Blanchet S: Beyond DNA: integrating inclusive inheritance into an extended theory of evolution. Nat Rev Genet. 2011, 12: 475-486. 10.1038/nrg3028.
PubMed CAS Article Google Scholar
49.
Rando OJ, Verstrepen KJ: Timescales of Genetic and Epigenetic inheritance. Cell. 2007, 128: 655-668. 10.1016/j.cell.2007.01.023.
PubMed CAS Article Google Scholar
50.
Bindra RS, Glazer PM: Genetic Instability and the tumour microenviroment towards the concept of microenviroment-induce mutagenesis. Mutation Res. 2005, 569: 75-85. 10.1016/j.mrfmmm.2004.03.013.
PubMed CAS Article Google Scholar
51.
Huang LE, Bindra RS, Glazer PM, Harris AL: Hypoxia-induced genetic instability - a calculated mechanism underlying tumour progression. J Mol Med. 2007, 85: 139-148. 10.1007/s00109-006-0133-6.
PubMed CAS Article Google Scholar
52.
Handel AE, Ramagopalan SV: Is Larmarckian evolution relevant to medicine?. BMC Med Genet. 2010, 11: 73-[[http://dx.doi.org/10.1186/1471-2350-11-73 ]]
PubMed PubMed Central Article Google Scholar
53.
Uhr J, Marches R: Dormancy in a Model of Murine B Cell Lymphoma. Seminars in Cancer Biol. 2001, 11: 277-283. 10.1006/scbi.2001.0383.
CAS Article Google Scholar

Download references

Acknowledgements

We thank very much the three referees, whose suggestions helped us to greatly improve this work.

MAJC gratefully acknowledges the support of the ERC Advanced Investigator Grant 227619, ’M5CGS - From Mutations to Metastases: Multiscale Mathematical Modelling of Cancer Growth and Spread’.

The AD’s work has been done in the framework of the Integrated Project “P-medicine - from data sharing and integration via VPH models to personalized medicine” (project ID: 270089), which is partially funded by the European Commission under the 7th framework program.

Author information

Affiliations

Division of Mathematics, University of Dundee, Dundee, Scotland, UK
Mohannad Al-Tameemi & Mark Chaplain
Department of Experimental Oncology, European Institute of Oncology, Via Ripamonti 435, Milano, I-20141, Italy
Alberto d’Onofrio

Authors

Mohannad Al-Tameemi
View author publications
You can also search for this author in PubMed Google Scholar
Mark Chaplain
View author publications
You can also search for this author in PubMed Google Scholar
Alberto d’Onofrio
View author publications
You can also search for this author in PubMed Google Scholar

Corresponding authors

Correspondence to Mark Chaplain or Alberto d’Onofrio.

Additional information

Competing interest

The authors declare that they have no competing interests.

Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

Authors’ original file for figure 1

Authors’ original file for figure 2

Authors’ original file for figure 3

Authors’ original file for figure 4

Authors’ original file for figure 5

Authors’ original file for figure 6

Authors’ original file for figure 7

Authors’ original file for figure 8

Authors’ original file for figure 9

Authors’ original file for figure 10

Authors’ original file for figure 11

Authors’ original file for figure 12

Authors’ original file for figure 13

Authors’ original file for figure 14

Authors’ original file for figure 15

Authors’ original file for figure 16

Authors’ original file for figure 17

Authors’ original file for figure 18

Authors’ original file for figure 19

Authors’ original file for figure 20

Rights and permissions

Reprints and Permissions

About this article

Cite this article

Al-Tameemi, M., Chaplain, M. & d’Onofrio, A. Evasion of tumours from the control of the immune system: consequences of brief encounters. Biol Direct 7, 31 (2012). https://doi.org/10.1186/1745-6150-7-31

Download citation

Received: 06 February 2012
Accepted: 26 July 2012
Published: 25 September 2012
DOI: https://doi.org/10.1186/1745-6150-7-31

Keywords

Tumour growth
Immune response
Cytotoxic T-lymphocytes
Immuno-evasion
Mathematical models
Chemotaxis
Diffusion
Immuno-editing

Evasion of tumours from the control of the immune system: consequences of brief encounters

Abstract

Background

Results

Conclusions

Reviewers

Background

Methods

The Mathematical Model

Spatiotemporal Dynamics of the Tumour Cells

Spatiotemporal Dynamics of the CTLs

Chemoattractant

Chemorepellent

Tumour cell-CTL Complexes

Modeling the transition rates and probabilities

Boundary and initial conditions

Results and Discussion

Spatially Homogenous Case

Spatiotemporal Model

Conclusions

Appendix

Non-dimensionalization

Numerical values of the parameters

Reviewers’ comments

Reviewer #1: Prof. G. Bocharov (nominated by Dr. V. Kuznetsov, member of the Editorial Board of Biology Direct) (Institute of Numerical Mathematics, Russian Academy of Sciences, Moscow, Russian Federation)

Reviewer #2: Prof. M. Kimmel (Rice University, Houston, USA)

Reviewer #3: Prof. A. Marciniak-Czochra (Heidelberg University)

The following were comments on the revised manuscript (other comments on missprints are omitted, and were implemented):

Author’s contributions

References

Acknowledgements

Author information

Affiliations

Corresponding authors

Additional information

Competing interest

Authors’ original submitted files for images

Rights and permissions

About this article

Cite this article

Keywords

Biology Direct

Contact us