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

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.

Moreover $k_0^{+}=1.3\times 10^{-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^{+}=\text{constant}=1.3\times 10^{-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^{+}=\text{constant}=1.3\times 10^{-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^{+}=\text{constant}=1.3\times 10^{-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 _i are 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:

$\text{MAX}\left(T_0\right)<\text{MAX}(T_1+\cdots +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^{+}=\text{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\left(t\right)=\underset{0}{\overset{1}{\int }}{T}_i(x,t)\mathrm{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 _i is 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^{+}=\text{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 .

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 _i is 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.

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.

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.