## Abstract

Using a hybrid cellular automaton, we investigate the transient and asymptotic dynamics of the cell-mediated immune response to tumour growth. We analyse the correspondence between this dynamics and the three phases of the theory of immunoedition: elimination, equilibrium and escape. Our results demonstrate that the immune system can keep a tumour dormant for long periods of time, but that this dormancy is based on a frail equilibrium between the mechanisms that spur the immune response and the growth of the tumour. Thus, we question the capacity of the cell-mediated immune response to sustain long periods of dormancy, as those appearing in recurrent disease. We suggest that its role might be rather to synergize with other types of tumour dormancy.

This article is part of the themed issue ‘Mathematical methods in medicine: neuroscience, cardiology and pathology’.

## 1. Introduction

Paul Ehrlich suggested in 1909 that the immune system might protect an organism from the development of cancer [1]. Around 50 years later, this proposition was more formally reintroduced by Macfarlane Burnet [2,3] and, later on, by Lewis Thomas [4]. After suffering major setbacks [5,6], the immunosurveillance theory gained renewed impetus close to 20 years ago, thanks to several experimental studies with genetically altered mice [7,8]. Currently, the immunosurveillance of tumours is more properly referred to as cancer immunoediting. Given the genetic heterogeneity of tumours, this control system coevolves with them and seems to act as a natural selective force, editing its phenotype by selecting those cells that are unresponsive to immune detection.

Cancer immunoediting can be described by three phases: elimination, equilibrium and escape. The first of these three Es [9] corresponds to what has traditionally been termed immunosurveillance [10], and involves the innate and the adaptive immune responses. During this phase, the immune system keeps in check a tumour cell population, successfully recognizing and destroying the majority of its cells. However, some residual tumour cells might remain unnoticed and asymptomatic for a long period of time, which can range from 5 years to more than 20 years. This period of time defines a second stage, in which a small cell population is kept at equilibrium. Finally, the phase of escape is led by some tumour cells that might present *a priori*, or have acquired along their evolutionary process, a non-immunogenic phenotype.

The mechanisms through which a tumour can be maintained at low cell numbers (i.e. dormant) are diverse. In a first approach, cancer dormancy can be generally classified in two categories: tumour mass dormancy and cellular dormancy [11]. In the former case, the equilibrium of a tumour is the result of a balance between cell growth and cell death. In the latter, the cells arrest and survive in a quiescent state until more benevolent conditions are provided by their environment. The occurrence of tumour mass dormancy is commonly associated with two different mechanisms [12]. The first is angiogenic dormancy, which occurs when the cells are unable to induce angiogenesis, and therefore to recruit oxygen and other nutrients to their location. In this manner, the proliferation rate is counterweighted by elevated rates of apoptosis. The second mechanism is the immune system response. This response is very complex and involves many types of cells and molecules [13]. There is evidence that the cell-mediated immune response collaborates with the humoral immune response to promote the dormancy of tumours, and that CD8^{+} lymphocytes and interferon-*γ* play a transcendental function in its maintenance [14].

In this work, we use a hybrid cellular automaton (CA) to investigate the dormancy of a tumour mass, mediated by the cellular immune response. Even though an interesting study has been previously carried out in this context [15], this study includes new features, which we believe makes it more realistic, permitting a correlation between the results and the theory of immunoedition. Mainly, the time scale of the cytotoxic cell action (about an hour) differs from the time scale of tumour cell proliferation (about a day). Secondly, our CA includes a new parameter that allows us to represent immunosuppressed environments. The exploration of different immunological scenarios enables the discussion of a possible dynamical origin of tumour dormancy and the sneaking through of tumours, as originally proposed by Kuznetsov *et al.* [16].

## 2. The model

We consider a model of limited nutrient growth of an immunogenic tumour, consisting of a hybrid CA. The model is very similar to another model presented by Mallet & De Pillis [15]. For a schematic representation of the model, see figure 1. These models are an extension of a previous CA model designed to study the effects of competition for nutrients and growth factors in avascular tumours [17]. They are hybrid because the tissue cells are treated discretely, allowing them to occupy diverse grid points in a particular two-dimensional spatial domain *Ω*, and evolve according to probabilistic and direct rules. On the other hand, the diffusion of nutrients required for growth and survival (such as glucose, oxygen or other types of nutrients) from the vessels, which are placed on the boundary of the spatial domain ∂*Ω*, is represented through linear reaction–diffusion equations, which are continuous and deterministic.

We utilize a quadrilateral spatial domain *Ω*=[0,*L*]×[0,*L*], which is partitioned into a regular grid with a resolution of *n*×*n* pixels, *n* being equal to 300 in all our simulations. Each of the grid points is occupied by one or more cells, that can execute several actions. This work includes four types of cells: healthy cells , tumour cells , immune effector cells and dead (necrotic and apoptotic) cells . Unlike previous studies, we do not make a distinction between the innate and the adaptive immune responses. For simplicity, we gather the natural killer cells and the T-cells in the same compartment, and simply refer to them as cytotoxic cells.

The role of the healthy cells is reduced to passive competitors for nutrients that allow the tumour cells to freely divide or migrate. The dead cells play no significant role in the model, and they can be replaced by the tumour and the immune cells, just as if they were phagocytized by wandering macrophages.

At each CA iteration, the tumour cells can carry out different actions attending to certain probabilistic rules. These rules depend on the nutrient concentration per tumour cell at a grid point and some specific parameters. Each of these parameters *θ*_{a} represents the intrinsic capacity of the tumour cells to carry out a particular action *a*. The tumour cells can divide *θ*_{div}, move *θ*_{mig} or die *θ*_{nec}. Attending to morphology, diverse types of tumours can be generated, depending on the nutrient competition parameters among tumour cells *α*,*λ*_{N}. Here, we consider parameter values that, in the absence of an immune response, generate rather spherical tumours [18].

Concerning the immune cells, we assume a constant input of cytotoxic cells into the region. When these cells are in contact with at least one tumour cell, they attempt to lyse it. The fate of this tumour cell depends on the intrinsic cytotoxic ability of the immune cells, which is represented by the parameter *θ*_{lys}. We outline the importance of this parameter, which allows us to represent immunosuppressed tumour microenvironments. Following previous studies [16,19], we assume that when a cytotoxic cell interacts with a tumour cell, several cytokines are secreted by the immune cells, which induce the recruitment of other immune cells to the domain. We note that the constant input of immune cells can, to some extent, also be regarded as a mechanism of recruitment. Commonly, cytokines are secreted by other types of immune cells, as for example T-helper cells, which are not explicitly modelled here. When immune cells are not in direct contact with a tumour cell, they can either move or become inactivated *θ*_{inc}. As in previous studies, we assume that a single immune cell cannot lyse more than three times, leaving the region of interest when this occurs [15]. The precise probabilistic laws are described in the following.

### (a) Diffusion of nutrients

Two types of nutrients are utilized in this model, making a distinction between those which are specific for cell division *N*(*x*,*y*,*t*) and others *M*(*x*,*y*,*t*) required for cell survival. The partial differential equations for the diffusion of nutrients are
2.1and
2.2where *T*(*x*,*y*,*t*), *H*(*x*,*y*,*t*) and *E*(*x*,*y*,*t*) are functions representing the number of tumour, healthy and immune cells at time *t* and position (*x*,*y*). For simplicity, we assume that both types of nutrients have the same diffusion coefficient *D*_{N}=*D*_{M}=*D*. Following [15], we consider that the competition parameters are equal *k*_{2}=*k*_{3}=*k*_{5}=*k*_{6}=*k*, except for the tumour cells, which compete more aggressively. We set *k*_{1}=*λ*_{N}*k* and *k*_{4}=*λ*_{M}*k*, with *λ*_{M} and *λ*_{N} greater than one. We exploit the difference in time scales for nutrient diffusion (minutes) and cell division (days), assuming that the solutions are stationary. On the vertical sides of the domain, where the vessels are placed, Dirichlet boundary value conditions are imposed. Therefore, we assign *N*(0,*y*)=*N*(*L*,*y*)=*N*_{0} and *M*(0,*y*)=*M*(*L*,*y*)=*M*_{0}. For simplicity, the horizontal upper and lower bounds of the domain obey periodic boundary conditions. Therefore, we impose *N*(*x*,0)=*N*(*x*,*L*) and *M*(*x*,0)=*M*(*x*,*L*).

The reaction–diffusion equations can be non-dimensionalized [17] by redefining the nutrients and the spatial and temporal coordinates in the form 2.3

Dropping the tildes and considering that the solutions are stationary, we obtain the equations
2.4and
2.5where *α*^{2}=*kL*^{2}/*Dn*^{2} is the dimensionless rate of consumption of nutrients by host and immune cells, while *λ*_{N}*α*^{2} is the rate of consumption of the nutrient *N* by the tumour cells. The boundary conditions now read *N*(0,*y*)=*N*(*L*,*y*)=1 and *M*(0,*y*)=*M*(*L*,*y*)=1. From a physical point of view, these elliptic partial differential equations represent the ‘scattering’ of nutrients from the boundary of a tissue and their diffusion in it. The role of the cells at a particular point in space is to act as a ‘potential barrier’, consuming nutrients and, therefore, attenuating their concentration at such position. The size of such ‘barrier’ varies in space, depending on the number of cells at each position and the rate at which they consume nutrients. As *λ*_{N}≥1 and *λ*_{M}≥1, tumour cells compete equally or more aggressively for nutrients than healthy cells.

### (b) Cellular automata rules

Now we describe the CA rules for the tumour and the cytotoxic cells. They are very similar to those used in [15], and are also described in [18]. Any difference will be explicitly remarked. In what follows, and are the tumour and the immune cells at position , while and are the concentration of nutrients in non-dimensional variables at position . We recall that the parameters *θ*_{a} represent the intrinsic capacity of a cell to carry out a particular action *a*.

#### (i) Tumour cell rules

Every 24 CA iterations the tumour cells are randomly selected one by one, and a die is rolled to choose whether each of these cell divides (action 1), migrates (action 2) or dies (action 3). When one of these three actions is selected, it might occur or not, depending on a probability distribution which depends on the nutrient concentration. The probability rule for each of these three actions is as follows.

(

*a*_{1}) A tumour cell divides with probability 2.6This probability is compared with the probability of a randomly generated number using a normal distribution and the same standard deviation. If the former is greater than the latter, division takes place. The higher the value of

*θ*_{div}, the more metabolic requirements for a cell to proliferate. When a cell at position divides, if there are neighbouring CA elements that are not currently occupied by tumour cells, we randomly select one and place there the newborn cell, thus making and or , where is the function representing the necrotic cells at position . However, if all the neighbouring elements are occupied, we let the cells pile up, making .(

*a*_{2}) A tumour cell migrates with probability 2.7If

*P*_{mig}is greater than the probability of a randomly generated number, migration proceeds, otherwise it does not. The higher the value of*θ*_{mig}, the more metabolic requirements for a cell to move, unless there are too many tumour cells. When a cell at position moves, if there are neighbouring CA elements that are not currently occupied by tumour cells, we randomly select one at and place the cell there. If there is more than one cell in the original position, the moving cell simply replaces the healthy or the necrotic cell, thus making the transformation , and or . On the other hand, if there is only one tumour cell at , then it interchanges its position with the healthy or necrotic cell at . If all the neighbouring elements are occupied, we displace the cell to a randomly selected neighbouring element.(

*a*_{3}) A tumour cell dies with probability 2.8If

*P*_{nec}is higher than the probability of a randomly generated number, necrosis proceeds, otherwise it does not. The higher the value of*θ*_{nec}, the greater the probability for a cell to die. When a cell at position dies, we make . If this is the only cell at , then .

#### (ii) Immune cell rules

At each CA iteration, there is a constant input of immune cells. These cells are placed at random in the domain, at points that are not occupied by tumour cells. Every such grid point is examined and, if a probabilistic condition holds, the healthy or dead cells that might occupy it are replaced with an immune cell (action 1). Then, the immune cells are randomly selected one by one. When the immune cell has one or more tumour cells as first neighbours, it carries out an attempt to lyse (action 2). If the tumour cell is destroyed by the immune cell, the first neighbouring cells are flagged for recruitment (action 3). Those effector cells whose immediate neighbourhood is not occupied by tumour cells either migrate or become inactivated. To decide which of these two processes is carried out, a coin is flipped. If the output is migration, it occurs for sure. In the opposite case, inactivation might result (action 4). The probability rules of these actions are as follows.

(

*a*_{1}) An immune cell is placed in the background with probability 2.9where*f*is a number between 0 and 1, that represents the intensity of the input of immune cells into the tissue. If*P*_{bkg}is greater than a randomly generated number between zero and one, then an immune cell appears in the corresponding grid point.(

*a*_{2}) An immune cell lyses a tumour cell with probability 2.10where*η*_{n}indicates summation up to the*n*th nearest neighbours. If*P*_{lys}is higher than the probability of a randomly generated number, then the selected tumour cell dies. Therefore, , and the immune cell counter decreases by a unit. If the counter (starting with a value of three) reaches a value of zero, the immune cell is inactivated and replaced by a healthy cell. The smaller the value of*θ*_{lys}, the greater the probability for an effector cell to lyse a tumour cell. This parameter was not present in [15] and is introduced here to model the intrinsic cytotoxicity of T cells.(

*a*_{3}) When an immune cell destroys a tumour cell, each CA element around it without tumour cells is explored, and a new immune cell is recruited with probability 2.11If

*P*_{rec}is higher than the probability of a randomly generated number, recruitment proceeds. The higher the value of*θ*_{rec}, the fewer the surrounding tumour cells that are required for T cell recruitment to succeed. When a cell is recruited at position , we make or , and .(

*a*_{4}) The inactivation of an immune cell occurs with probability 2.12If*P*_{inc}is higher than the probability of a randomly generated number, inactivation proceeds. The smaller the value of*θ*_{inc}, the fewer the surrounding tumour cells that are required for a T cell to become inactivated. When a cell disappears from position , we simply make and .

### (c) The algorithm

The algorithm starts with a domain full of healthy cells, except for a single tumour cell placed at the centre of the domain. Firstly, we let the tumour grow until it is detected by the immune system, when it has reached some specific size *T*_{det}. During this period of growth, each CA step corresponds to one day. Each iteration begins with the integration of the reaction–diffusion equations, using a finite-difference scheme and a successive overrelaxation method. Then, all the tumour cells are randomly selected with equal probability, and the CA rules are applied. As in previous studies [17], every time an action takes place, the reaction–diffusion equations are locally solved in a neighbourhood with size 20×20 grid points. When the time of detection is reached, the immune cells start to evolve. Now the CA step corresponds to 1 h, and during the next 23 steps, only the immune cells are computed. First, the background source of immune cells is executed. Then, the reaction–diffusion equations are solved and all the immune cells are randomly selected. For each immune cell, after applying the CA rules, the nutrients are computed in a local region, in exactly the same manner as before. Every 23 iterations, the tumour cells are checked and the tumour cell rules applied as previously described, before immune detection. The algorithm stops when a maximum number of steps after the elapse of the immune response has been reached, or when a tumour cell is at a distance of two grid points from its boundary.

## 3. Simulations and results

We study the evolution of the tumour and the immune system for three different scenarios. The first scenario is used as a reference, and it is characterized by high levels of immune cell recruitment and negligible necrosis due to the scarcity of nutrients in the core of the tumour masses. In the second scenario, the recruitment levels are reduced, while the necrosis of tumour cells is enhanced in the third. Unless specified, the remaining parameters are all the same in every case. Beginning with one tumour cell, the tumours grow up to *T*_{det}=5×10^{3} cells, and at this moment the immune response triggers. In order to elucidate the effects of tumour immunogenicity, we devise what shall be called a *transient bifurcation diagram*. Given a dynamical system, a bifurcation diagram is a plot of the asymptotic values of a particular variable against a set of values of some relevant parameter. However, in many situations there might exist very long transients before the asymptotic state is established. These transients are of great importance in our context, because tumours may exhibit long periods of latency before the development of recurrence. Therefore, we compute the number of tumour cells for the last 100 iterations of a trajectory comprising 24 000 iterations of the CA from immune detection. Then, these 100 points are represented on the vertical axis for different values of the parameter *θ*_{lys}, which lies on the horizontal axis. If we assign to each of these iterations a time of one hour, we are registering the size of the tumour for approximately the last 4 days of a period of 33 months from immune detection. We recall that the parameter *θ*_{lys} codes the intrinsic ability of the immune cells to recognize and lyse their adversaries. Higher values of this parameter correspond to more immunodeficient environments.

### (a) Reference scenario

The set of parameters for this scenario is chosen similar to previous studies, in which it has been demonstrated that they generate reasonable tumour dynamics [17,18]. The specific values are *θ*_{div}=0.3, *θ*_{nec}=0.05, , *θ*_{rec}=1.0, *θ*_{inc}=0.1, *λ*_{M}=10, *λ*_{N}=25 and *α*=2/*n*. Regarding the natural flow of immune cells into the tissue, two situations are inspected for each scenario. The first corresponds to a high input of immune cells into the tumour area. In this case a value *f*=0.10 is set, which means that approximately 10% of the background is occupied by immune cells, if there are not too many immune cells piled up. The other has a lower input of 5%, thus *f*=0.05. In the absence of immune response, the tumours display a rather spherical shape. As we can see from the transition bifurcation diagrams shown in figure 2, three different regions are clearly distinguished. In the first region, when the immune system is effective, the tumours are completely eliminated. The second is related to an equilibrium phase, for which tumours spend very long transients oscillating at low cell numbers. Finally, tumours with increasing size, eventually leaving the domain through the vessels, appear in the third region. Thus, here we see how immunogenicity affects the fate of tumours, in accordance with the theory of immunoedition.

To give insight into the second and the third regions, time series have been computed (figures 3 and 4), until the tumour escapes. Initially, the tumours grow in the absence of immune response, and then the immune cells start to reduce them or, in the worst case, delay their growth. Depending on how effective the immune cells are, longer or shorter transients follow this reduction phase. The asymptotic dynamics is always the same: if the tumours are not totally eliminated by an efficient immune system, they eventually escape from the region. These two attractors are reminiscent of those appearing in [16]. As shown in figure 3*a*, the length of the transients in the second region, which are of around 12 years, clearly indicate a phase of prolonged tumour mass dormancy. During the period of dormancy the immune system keeps the tumour at low cell numbers and randomly displaces its disconnected pieces until one of them reaches the vessels. In the third region, transients are found again, but they are shorter (less than 4 years) and the tumours at escape have bigger sizes. As predicted by Kuznetsov *et al.* [16], the duration of the transients is stochastic. This randomness is evident from the transient bifurcation diagrams, since after 33 months of tumour-immune struggle, some tumours have escaped and some others have not, disregarding how immunogenic they are. When the immune system barely responds to the tumour, we see very big tumours occupying the domain and escaping rapidly, as depicted in figure 3*e*.

Interestingly, the equilibrium region shrinks as the normal input of cells into the tissue reduces from *f*=0.1 to *f*=0.05. As is shown in figure 4, the oscillations during the equilibrium phase are more pronounced. This makes the equilibrium more unstable and suggests that having cells scattered all over the domain is important for the maintenance of dormancy. Probably, the reason is that these spread immune cells keep the tumour at a small size, not allowing its overgrowth in any specific direction.

We have also explored the importance of the tumour size at detection by reducing this size to 5×10^{2} cells. The results are depicted in figure 5 and resemble very much those shown in figure 2. There is no hint of a sneaking through mechanism in our model. According to the definition given by Gatenby *et al.* [20], sneaking through is the preferential take of tumours after small size inocula to a similar degree to that seen with large size inocula, compared with the rejection of medium sized inocula. More clearly put, small and big tumours escape immune surveillance, while intermediate ones do not. Such phenomenon has not been observed in the present case for other values of the tumour size at detection. However, we do not discard it, because motility of tumour cells has not been included in this first investigation, and might be crucial for these cells to escape.

Finally, even though the tumours here inspected are genetically homogeneous and no evolutionary process is really taking place in our model, the transient bifurcation diagrams insinuate how the sculpting of the phenotype occurs, moving from the first region to the second, and then to the third. In fact, a similar CA can be used to explore the impact of heterogeneity and how the process of immunoedition takes place. It suffices to consider that the immune cell intrinsic cytotoxicity, represented by the parameter *θ*_{lys}, depends on the tumour cell.

### (b) Low recruitment scenario

We now evaluate the impact of the recruitment of immune cells to the domain of the tumour. For this purpose, we reduce the value of *θ*_{rec} from 1 to 0.35. Our interest in this parameter is due to the fact that, on many occasions, the recruitment of cells to the site of the tumour might be very complicated. The recruitment of immune cells is a very complex process, at least from a physical point of view. The extravasation of leucocytes requires an initial contact between these cells and the endothelial cells, which depends on adhesion molecules. After adhesion to the walls of the vessels, the immune cells traverse them through diapedesis, which again relies on several cytokines. Finally, chemokines bias their random walks to the tumour location [21]. Thus, we expect this parameter to exhibit great fluctuations, depending on the tissue location and other factors, as for instance the degree of inflammation.

The effects of decreasing the recruitment parameter are shown in figure 6. As expected, the elimination region shrinks, while the escape region widens. A dramatic reduction of the dormancy window is observed in both plots. When *f*=0.1, the window still exists, but for *f*=0.05 it has disappeared. These results suggest that a relatively tight balance between lysis and growth is required to maintain the dynamical equilibrium of the tumour.

Note that, as previously proposed, the equilibrium of the tumour implies that reduction must occur in an isotropic manner. If a region of the tumour grows over the immune system capacity, then a soon overgrowth and a consequent escape would be expected. In this model, this relies on a positive feedback mechanism between the natural input of immune cells and their recruitment. The more cells there are spread in the domain, the more chances for an immune cell to lethally hit a tumour cell. When this occurs, recruitment proceeds, favouring the local aggregation of immune cells at this site of the tumour and giving rise to satellites [15]. This isotropy can be appreciated in the equilibrium represented in figures 3*b* and 4*b*, as opposed to those situations that lie in the third region, represented in figures 3*d*,*f* and 4*d*.

### (c) High necrosis scenario

Solid tumours exhibit sometimes necrotic cores due to the scarcity of nutrients. Other chemical species can be represented with the present model (e.g. growth factors) and, if desired, necrosis can be regarded as apoptosis, at least to some extent. Therefore, we now inspect the effects of cell death in the model. To this end, we increase the value of *θ*_{nec} from almost zero to 0.5. Obviously, the increase of necrosis facilitates the labour of the immune system. As shown in figure 7, the elimination region enlarges substantially, compared with the reference case. Also in the equilibrium region, lower tumour cell numbers are seen before the escape of the tumour. More importantly, the equilibrium window, which has been associated with large periods of tumour mass latency, is practically imperceptible for *f*=1.0 and has completely vanished for *f*=0.05. We have again computed time series, showing that transients occur in the equilibrium region, sometimes as long as those appearing before in the equilibrium, but generally shorter (figure 8). In fact, the equilibrium window and the escape zone drawn in figure 7*a* overlap. It seems that the equilibrium region appearing in the reference scenario has been swept under the elimination region. Once more, the results confirm the requisite of a relatively delicate balance between the mechanisms that maintain the cytotoxic destruction of the immune system and the growth of the tumour, in order to keep it at low cell numbers for long periods of time.

## 4. Conclusion

We have studied the transient and asymptotic dynamics of a CA model for tumour–immune interactions. We have shown that, depending on the immunogenicity of the tumour, the model furnishes three main types of dynamics, which are in close relationship with the three phases of the theory of immunoedition. Importantly, we have shown that a dynamical equilibrium between the tumour can occur for long periods of time, as proposed by Kuznetsov *et al.* [16]. However, after inspection of the parameter space, our model suggests that this equilibrium is quite fragile, because it is based on an adjusted balance between the mechanisms that stimulate the immune response and tumour cell proliferation. It is worth asking if this also occurs in the model presented by these authors [16]. We must bear in mind that such a model is very sensitive to a parameter, there called *μ*, which is related to the rate at which tumour cells are lysed. In their model, tumour cell lysis is represented by a simple Lotka–Volterra competition term. In other studies [18,19,22,23], it has been demonstrated that the rate at which a tumour is lysed is represented better by a more sophisticated function, that depends on the geometrical properties of the tumour and its immunogenicity. Moreover, their work did not assess the importance of the parameter that models the recruitment, there represented by the Greek letter *ρ*. Although their model is tested against experimental data obtained from a BCL_{1} lymphoma in chimeric mice, these tumours develop in the spleen of the mice. The spleen is a secondary lymphoid organ through which T-cells are permanently trafficking, and the process of recruitment to other types of tissues might be more arduous. Interestingly, a reduction of the value of the parameter *ρ* from 1.131 to 0.630, which in their model is related to the rate at which T-cells are recruited, produces a saddle-node bifurcation through which the dormant state disappears. Just as in this work, considerable levels of recruitment are required to sustain dormancy. Furthermore, the infiltration of the immune cells into the tumour mass has been neglected in this work. Presumably, this effect would make the equilibrium even more delicate. Nevertheless, both models clearly demonstrate that a state of tumour mass dormancy mediated by the immune system can occur. It is the length of this dormant period that it is being questioned here. Thus, we conclude that, even though tumour mass dormancy can result from the cell-mediated immune response to tumour growth, long periods of dormancy, as commonly found in recurrent metastatic tumours [11,12], are not likely to arise by this single mechanism. It is therefore pertinent to ask ourselves if the role of the cell-mediated immune response in the promotion of the dormancy of a tumour mass is rather to synergize with other types of more efficient mechanisms, as for example cellular dormancy.

## Authors' contributions

All the authors have contributed equally to this work.

## Competing interests

We declare we have no competing interests.

## Funding

This work was supported by the Spanish Ministry of Economy and Competitiveness under project no. FIS2013-40653-P and by the Spanish State Research Agency (AEI) and the European Regional Development Fund (FEDER) under project no. FIS2016-76883-P. M.A.F.S. acknowledges the jointly sponsored financial support by the Fulbright Programme and the Spanish Ministry of Education (programme no. FMECD-ST-2016).

## Footnotes

One contribution of 14 to a theme issue ‘Mathematical methods in medicine: neuroscience, cardiology and pathology’.

- Accepted November 7, 2016.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.