The role of interventions in the cancer evolution–an evolutionary games approach

1.   Introduction

The objective of this study is to develop and apply a new mathematical model to address the important questions for basic cancer research and practice of oncology: What is the role of different kinds of intervention (exposures or treatments) in the evolution of cancer? Can proliferation acceleration make carcinogenesis more likely? Can cell killing create conditions that enable the evolution of heterogeneity and therefore increase the viability of the cancer cell population? An essential part of the problem we address is whether the genetic heterogeneity in cancer cells exists ab initio (before diagnosis and intervention), or does it evolve as a result of the intervention. We believe, basing on evidence, that the former is the case, as stated in the following hypothesis: positive (stimulating) or negative (cell-killing) treatment modifies pre-cancer and tumor cell heterogeneity through the promotion of selection of preexisting cell clones, their phenotypic plasticity and self-organization. We are far from believing that using the new model we prepared will prove completely this hypothesis. Nevertheless, we hope that the new approach based on the theory of evolutionary games gives an efficient tool for investigation of different aspects related to the hypothesis. The evolutionary game theory (EGT) combines machinery of the theory of games with present knowledge of population biology and evolution ^[28,30]. EGT differs from standard game theory in deviating from the rational approach of the competitive players, in considering treatment of strategies as phenotypes of individuals acquired through the evolution. Moreover, the players are members of a population consisting of individuals with different phenotypes (strategies), who can cooperate or compete for resources. As a result of various adaptations to the environment and following games in time (generations) the population can tend to stabilize its structure at the same time gaining stable monomorphism or polymorphism of population's phenotypes. Such state is called evolutionary stable. Whereas evolutionary stable strategy (ESS) is defined as a phenotype that, if adopted by the vast majority of a population, will not be displaced by any other phenotype ^[29]. However, opposite situation is very probable to occur.

Classical models of cancer development assume that mutations, which promote the growth of cancer cells affect only the units in which they occur ^[13,23,35] but recent studies point out that tumor cells are able to adopt various genetic strategies which may influence the rate of their own development. Moreover the mutation can also affect neighboring cells ^[37], and combined with cooperative behaviours, competition for resources such as space, oxygen and nutrition, occurs between different subspecies within the same tumor ^[1,26] in which internal communication between tumor cells, and between tumor and normal cells, their competition for resources, hierarchical subordination, and collaboration, play an important role in cancer development and differentiation or disease transmission and reaction to stress factors including therapy ^[4,15]. The evolutionary game is performed between cells with different phenotypes (both normal and cancer cells). The main aim of these game theoretic models is to study the possibility of coexistence or even domination of newly formed tumor cells, which have acquired new strategies (phenotypes) by mutations. To our knowledge, Tomlinson and Bodmer ^[37] first proposed such a model describing inter-cellular interactions including avoidance of apoptosis and production of angiogenic factors. The models that followed their research described phenomena such as: production of the cytotoxic substances ^[36], production of growth factors ^[3], invasion and metastasis ^[25], tumor-environment interactions ^[14], radiation bystander effect ^[21], resistance to chemotherapy and p53 vaccine ^[7], interaction between osteoclasts and osteoblasts ^[11], tumor-stroma interaction ^[15], interaction between different tumors ^[5] and others (see ^[4,10,17,32] for survey). Usually, the analysis is limited to two and three phenotypes. Our previous paper ^[20] makes for the exception in which interaction between four different phenotypes of cells are illustrated using three-dimensional simplexes and time courses. The work of Basanta et al.^[8] is the only other paper in which interaction of four phenotypes was discussed. The authors, however, were interested in changes in subpopulations of chosen phenotypes with respect to changes of cost parameters rather than in studying equilibrium between all phenotypes and dynamics of their evolution. The reason for that is related to difficulties with analysis of simplexes of dimension higher than three and complexities in dynamics represented by nonlinear replicator equations increasing with their dimension. It is important to notice that dimension of replicator dynamics equations in the case of three phenotypes is equal to two which means that complex dynamical behaviours typical for nonlinear dynamics should be absent. In our opinion, it is one of the major disadvantages of the small number of considered strategies. Using EGT, we are able to predict whether the population has a tendency to become heterogeneous or rather only one phenotype will survive dominating the overall population. Moreover, using the so-called replicator dynamics ^[16] we may track the rate of population in time starting from the initial state to the equilibrium defined by ESS. Nevertheless, because of assumptions that perfect mixing occurs it gives only mean field results without possibility to take the effect of local arrangements on intercellular interactions into account. This drawback could be overcome by the methods of spatial evolutionary games (SEGT) based on cellular automata which are the next step in the discovery of new behaviors among cells and give different results than mean-field (averaged, "well mixed", non-spatial and giving values for the population without spatial information ^[22]) models. Nowadays, spatial games are quickly gaining popularity. Nevertheless, it should be remembered that the use of cellular automata ^[27] in conjunction with the classical theory of games originated spatial games. The line of reasoning presented by Bach et al. ^[2], where spatial tools used in modeling of carcinogenesis is best suited to our expectations. Similarly to non-spatial games, the spatial ones are also iterated. We proceed transient generations according to the following steps ^[2]: payoff updating - a sum of local fitness of neighboring cells, cell mortality - removing a certain number of players, reproduction by competition - defining which of the cells (or strategies) will be on an empty place. The game is played on the lattice forming torus, and all the competition results giving tie are settled randomly. The authors ^[2] present three ways of cell mortality: synchronous updating - all the cells die simultaneously, and they replacement depends on the strategy of their neighbors before dying, asynchronous updating - in every generation a single cell, chosen at random, dies and is replaced, semi-synchronous updating - probability of individual cellular mortality is equal to 0.1. In the last type 10% of players are deleted from the lattice in every generation. In our studies we are using mainly semi-synchronous updating, since this method allows for the biologically realistic situation. Furthermore simulations show that synchronous updating assumes a global controller of the system, while asynchronous updating implies that the vanishing of small cell clusters is impossible. The reproduction of the removed players (killed cells) is the next stage of the algorithm. Bach et al. have suggested two kinds of reproduction: deterministic - in competition for the empty place the winner is the strongest player (with highest local adaptation - sum of eight scores from cell-cell interaction), probabilistic - values of adaptation (or sum of values from payoff matrix) for each player are divided by total score in their neighborhood. Such local competition, with appropriate fitness and location, allows for strategies of cells with lower fitness, but in better location and locally superior in numbers to dominate in the population. In our opinion ^[34] deterministic reproduction describes a response to direct exposure of cells to an external stress factor while the probabilistic one is more adequate in modeling a bystander effect in cell populations. The main assumption of the spatial games presented is that each cell on the lattice is represented by a player following only one strategy. The local payoff for each player is the sum of payoffs due to interactions (according to the payoff matrix) with cells in the neighbourhood. We will refer to this approach as the standard, classical one or SEGT. Cells on the spatial lattice can also be considered as heterogeneous (instead of homogeneous), so that each particular player may contain mixed phenotypes. Spatial games of the type proposed in ^[19,20,33] are called mixed spatial evolutionary games (MSEG). It is important to mention the definition of the phenotype, which is the set of traits or characteristics of an organism ^[24]. Hence, in MSEG different degrees of playing a particular strategy are treated as different characteristics that define different phenotypes. It may happen that within the population, all of the players have diverse phenotypes (which probably better describes biological phenomena). For the sake of simplicity and following the way of reasoning from SEGT, those strategies and traits still correspond to the phenotypes and a general, collective point of view is defined as a players' phenotypic composition. In fact, the game is performed on a multidimensional lattice (dependent on the number of defined phenotypes in the model), where each layer represents a particular phenotype (as the frequency of occurrence) of the player. For the computation of the local adaptation, the sum of the payoffs between each phenotype (within two players) multiplied by their rate of occurrence is calculated first. The second step is summing of these values for each player in the neighbourhood. Another difference between SEGT and MSEG is that the tie (when payoffs are equal) for the former is settled randomly, while for the latter the average of phenotypic compositions is computed. The payoff matrix does not have to remain constant during the game; the parameters may change and may depend on other factors (e.g. changes in the environment or any other external influence). The model may be extended by yet another parameter $r$, which represents the amount and the availability of resources like territory, food, drugs etc. For mean-field games, $r$ is changed for the entire population at once and does not depend on the current consistency of the population. Spatial games give an opportunity to change $r$ for a group of cells or even for each cell separately. For the spatial application, we choose MSEG games since they provide more possibilities of defining different rules for the resources simulation. In this way, another lattice is needed, where a different game with different rules is performed. To explain the role and use of resources in the evolutionary game model we start with a classical hawk-dove model analyzed by mixed spatial evolutionary games (MSEG) supported with an additional spatial layer representing changes within the environment. Then we use the same approach to study the so-called angiogenic games ^[36], in the version containing varying external resources. The last model considered in this study is the previously mentioned four phenotype game ^[20] with two resources representing two different strategies of external interventions.

2.   Hawks and doves with resources

Hawk and Dove game is one of the first evolutionary models proposed by John Maynard Smith ^[30]. It includes two kinds of behaviour, fight or avoid, within a population of one species which is a symbolic representation of the ritual conflicts of two different strategies that arose in the process of evolution.

This game is an example of a two-player game, where each player can choose a strategy from a finite set of strategies. This game has two players, participants of the game and every player have its own decision set, which we call here a strategy. Every pair of players strategies will result in some game outcome for each player, which we call a payoff. Those are the values, written in the matrix form, that can be treated as a profit or cost of selecting a particular strategy. This values can also model, for example, Darwinian fitness. The payoff matrix (Table 1) contains two parameters: $v$ - the benefit that may be obtained during the contest and $c$ - the cost of escalation and it should be read horizontally, for a row-player way.

For non-zero sum game, in a bimatrix game where payoff matrices are respectively $A = \{a_{ij}\}$ and $B = \{b_{ij}\}$, a noncooperative equilibrium is defined by a pair of strategies $(a_{i^*j^*}, b_{i^*j^*})$ for which the following pair of inequalities is satisfied:

for all $i = 1, ..., m$ and $j = 1, ..., n$ (assuming that players are maximizing their payoffs). In the symmetric two-player game, like Hawk and Dove, $B = A^T$, so this condition can be simplified to $a_{i^*j^*} \geq a_{ij^*}$; $a_{j^*i^*} \geq a_{ji^*}$. There are cases, however, when pure Nash equilibrium strategies might not exist or there may exist many admissible equilibria. One may introduce, the so called, mixed strategies, where the game is played many times in the same conditions. The mixed strategy is in this case probability of playing pure strategies. For the Hawk and Dove game, the average payoff for Hawk will equal: $x(v-c)+(1-x)2v$ and for Dove: $v(1-x)$, if $x$ is the probability of playing Hawk.

Extending the notion of repeating game for the whole population, where we consider many players, that play the two-player game, we will reach the idea of evolutionary games. Here, the role of the individual player is negligible and the resulting payoff depends on the frequency of strategy played by the members of the population. This value characterizes the population and each individual interacts with other using some phenotype rather than choosing any strategy. For those games we can determine the evolutionarily stable strategy (ESS) $X$ that satisfies conditions that for any mutant type $p$ we have: $E(X, X) \geq E(p, X)$ (Nash equilibrium) and if $E(X, X) = E(p, X)$, then $E(X, p) > E(p, p)$. It means, that if in the population where strategy $X$ is presented by individuals, some mutant $p$ will occur, he will not be able to dominate the $X$ population, but the opposite situation is possible. In EGT the dynamics can be described by the replicator dynamics equation (RDE), in general form :

where: $e_i$ - is the individual strategy, $x$ - is the mean strategy, and E() - is the average payoff function.

In this particular case for the Hawk and Dove game we obtain:

For Hawk and Dove game stable polymorphism (coexistence between all phenotypes) is defined by the following ESS: $(v/c, 1-v/c)$ if $v<c$. This conditions could be found directly from the definition of ESS or from the Bishop-Canning (B-C) theorem ^[9]. The use of B-C theorem is not necessary but in the case of greater number of phenotypes it seems to be efficient method of finding ESS. One can easily check that this is an asymptotically stable equilibrium of the replicator dynamics equation (3). Generally, each ESS should have this property but reverse is not true (there may exist an asymptotically stable equilibria of RDE which are not ESS). Hence, to achieve a stable polymorphic result $v$ must be less than $c$, otherwise the population is dominated by Hawks. The results for this model are independent of the initial frequencies of occurrence (Figure 1).

Now we introduce a new parameter $r$ which represents additional resources available for players. We define some general rules:

● Hawks, being more agammaessive, tend to be better adjusted when $r$ is relatively small.

● Doves do better (or at least not worse than in a classical game) when resources are greater, providing opportunities, for instance, to "steal" food or to find it without interaction with another individual in the vast area.

● Even if $r$ is relatively large, hawks still bear the cost of escalation and doves still share the benefits. The exact value of $r$, however, may affect the intensity of escalation or distribution of resources.

The parameter $r$ is time varying, for example, it can be described by a piecewise linear function or as a step function (changing from minimum to maximum value or in the opposite direction). In this case the stationary solution can exist. However, if we assume an arbitrary function r(t) than the stationary solution does not exist, generally. The attractor may not exist or can be very complex. In particular case, where function r(t) is simple, the point attractor can be found analytically or by simulations.

Thus the payoffs for particular interactions are now functions dependent on $r$, which are linear and are presented in the form of payoff values for $r$ equal 0 and 1 (Table 3). Red font indicates the changes in comparison with original model (Table 1). The assumption of this model is that for $r$ = 0 the payoff matrix is the same as that in the basic model. With increasing $r$, more chances and better adjustment are gained by doves (observe red numbers in Table 3). Therefore doves do not have to share the benefits amongst their own population, and in the case of interaction with hawks it is still possible to gain some profits. Since the game is of a mean-field type, thus $r$ is changed for the entire population at once and does not depend on the current consistency of the population.

Using the RDE (3), we can formulate the dynamics of this new model including parameter $r$. In this case, the previous equation (4) will take the following form:

which for $r = 0$ gives (4).

Spatial games give an opportunity to change $r$ for a group of cells or even for each cell separately. All spatial simulations in this paper (for all presented here models) were calculated for a 2D square lattice witch borders (30 x 30 subjects) are neighbouring, that is why it has a shape of a torus. For the spatial application we choose MSEG games since they provide more possibilities of defining different rules for the resources. The standard assumption of the spatial games (called here SEGT) is that each cell on the lattice is represented by a player following only one strategy. The local payoff for each player is the sum of payoffs due to interactions (according to the payoff matrix) with cells in the neighbourhood. However, if we consider cells on the spatial lattice as heterogeneous, then each player may contain mixed phenotypes (called here MSEG - mixed (multilayer) spatial evolutionary games). In this case the strategy played by the cell is almost arbitrary depending on a number of unknown environmental conditions. The choice of a particular strategy may result in cell differentiation and escape to the population of differentiated cells. Hence, in MSEG different degrees of playing a particular strategy are treated as different characteristics that define different phenotypes. It may happen that within the population all of the players have diverse phenotypes (which probably better describes biological phenomena). In fact, the game is performed on a multidimensional lattice, where each layer represents a particular phenotype (as the frequency of occurrence) of the player. For the computation of the local adaptation, the sum of the payoffs between each phenotype (within two players) multiplied by their rate of occurrence is calculated first. The second step is the summing of these values for each player in the neighbourhood.

First we start with the mean-field model where parameter $r$ is changed piecewise linearly from 1 to 0. The aim is to show the possible results without focusing on the $r$ profile during the following generations. The set of parameters used for Figure 1 shows a population where the doves initially dominate, but then while resources decrease it impacts also the number of doves. We can also see some kind of dynamic delay : the frequency of occurrence of doves still rises a little after the resources have declined. When the benefits are higher than costs ($v$ = 12, $c$ = 6), the dynamics of changes depend on the initial frequencies of occurrence at $r$ = 1. Moreover, the frequencies are constant for that period (Figure 2). Still the final population contains only hawks.

For the MSEG model we introduce additional rules:

1. Individuals who are more hawks take an amount of resources from those who are less agammaessive, which amount results from the difference between hawks frequency of occurrence.

2. A dove shares its resources with those doves that have less resources according to the dove's frequency of occurrence.

3. The resources are limited to the (0, 1) interval.

An example describes this idea in a clearer way. The players are: $A - H = 80\%, D = 20\%, r = 0.7, B - H = 30\%, D = 70\%, r = 0.5$ Player A: the first rule requires that 50% of player B's resources shall be added to player A. Thus in the next iteration $r$ for player A is equal to 0.95 and for player B to 0.25. The second rule does not have an application here. If the difference between hawks is larger and for instance 0.4 of the resources should be transferred from player B to player A, then this value shall be limited to 0.3 since $r$ can be maximum 1 (hence, player B loses only 0.3).

Player B: the first rule does not have an application here, since it is dependent on the frequency of occurrence of hawks. Due to the fact that player A has got more resources, rule 2 can be applied. The difference between resources shall be multiplied by 70%, which gives 0.14 to be added to player B ($r$ = 0.39) and subtracted from player A ($r$ = 0.81). The model with an additional dimension (Table 3) will lead to increased adjustment for doves, especially for high values of resources. Therefore for $v$ = 3, $c$ = 9 and $r$ = 1 doves totally dominate in the population for the mean-field model (if $r$ does not change). Figure 3 shows the changes of particular phenotypes in the entire population for the model with and without resources. We can observe differences in both probabilistic and deterministic reproductions. The dynamics of changes within the probabilistic one oscillate more than in the basic (without resources) MSEG and hawks are a bit better adjusted (see upper row in Figure 3). Major differences are seen within the deterministic reproduction (which gives a similar result as the probabilistic one) where doves no longer dominate in the population. Introducing resources in hawk-dove model (Table 2) shall result in dove domination (Table 3-doves are better adjusted when $r>0$ comparing with classic model). Surprisingly, hawks have better adjustment than in the basic MSEG model. Due to the fact that starting from $r = 0$ the payoff matrix is the basic one and when $r = 1$ then doves shall dominate, the expected result is that doves shall be the dominant phenotype in the population.

For the set of parameters $v$ = 12 and $c$ = 6 (Figure 4) only probabilistic reproduction gives similar results to those for MSEG without resources.Deterministic reproduction shows some discrepancies in contrast to not enhanced MSEG. Moreover, interesting clusters of different cells are feasible what is shown in Figures 5 and 6.

3.   Angiogenic game with resources

One of the first models based on the evolutionary game theory describing results of interactions between tumor cells was proposed by Tomlinson and Bodmer ^[37] and described the paracrine production of growth factors including angiogenic promoters. The two phenotypes considered in this model represent cells producing growth factors (in paracrine fashion) and cells which do not produce them, and model parameters are related to costs of proangiogenic factor production (parameter $i$ in the payoff matrix, see Tables 4 and 5) and benefit resulting from these factors (parameter $j$ in the payoff matrix, in Tables 4 and 5). To reach stable dimorphism between the two phenotypes, the costs should be greater than the benefits, and the resulting frequencies depend on the ratio of these parameters. In the opposite case, the cells which do not produce growth factors dominate the population, and their phenotype is the only evolutionary stable strategy. A spatial version of this game ^[32] leads to similar effects. An extension of the angiogenic game was proposed by Bach et al. ^[3]. The authors assumed that to release the growth factors some synergy is needed and benefit from growth factors was obtained only when at least two out of three cells interactions in the same time represent the phenotype responsible for the production of proangiogenic factors. In this case, the stable dimorphism is possible only when the benefits are at least twice as great as the costs, but it also depends on the initial distribution of phenotypes in the population. A model with resources we propose takes an external intervention into account represented as treatment of the population using proangiogenic factors. Payoff tables for the standard angiogenic game and its version with resources are given by Table 4 and Table 5. In this case, we assume that resources represented by $r$ may change from 0 to 2. Results of simulations of MSEGs, as well as resources distribution, for different types of reproduction are presented on Figures 7 and 8.

The results qualitatively are similar to those for the hawk-dove game. An interesting, but quite intuitive finding is that additional external growth factors result in changes in final distribution and transient dynamics especially in the case of the deterministic reproduction. Production of angiogenic factors becomes less profitable than in the case when there is no external intervention. Moreover, the spatial distribution of the use of resources by the cells remarkably depends on the initial distribution of phenotypes on the lattice.

Similarly to the previous model we can write the replicator dynamics equations (Eq. 6), in the following form for the original angiogenic model as well as for the model with resources:

The introduction of resources in this model will change the baseline for the payoff matrix, which, in model without resources, was equal to 1 in our case. Nevertheless, since $r$ is varying in time it causes that the "new" baseline cost of meeting is also time-varying. Additional factor in the payoff, that comes from resources, added to each element of the payoff matrix will not change the dynamics of the system, that is why the replicator dynamics is the same for both models - with and without resources. The difference in the outcome of the model with resources is a result of non-homogeneous resources distribution. An interesting finding is that the clusters of phenotypes in the case of probabilistic reproduction are almost "reproduced" by clusters in resources (Figure 7-upper left (b), and Figure 8 left (b)). When condition for existence of point of equilibrium is fulfilled ($i<j$), then we can observe the coexistence of both populations as it is shown in Figure 9. In this case $i = 0.4$ and $j = 0.6$ and in the final MSEG lattices some shades of blue color can be observed, which is a result of mixing with a little amount of green ($A-$ phenotype).

4.   Four phenotype model with resources

To illustrate advantages of our approach to the analysis of increasing number of strategies let us consider the model which combines two classical models of Tomlinson ^[36,37], and one of them-the angiogenic game-was described in the previous section. The model contains four different strategies/phenotypes of cells:

● A cell produces a growth factor and the benefit impacts all the neighbors and the cell itself (phenotype A);

● A cell produces a cytotoxic substance against nearby cells (phenotype P);

● A cell is resistant to the cytotoxic substance (phenotype Q);

● Strategy which shall be considered as a baseline (no production of the cytotoxic substance, no resistance to it, no growth factor) (phenotype R).

In this case, the proposed pay-off table is given by Table 6 and its counterpart with resources by Table 7. We assume the existence of two different resources one of which represents an effect of growth factors ($r$), and the other killing effect ($c$).

In the first case, the external growth factor can have a positive (to recruit cells from the dormant or quiescent phase) or negative effect (as angiogenesis inhibitors, for example, Sunitibib in antiangiogenic therapy). The way they appeared in the payoff matrix results from their actions. Hence the similarity to the cytotoxic and growth effects of the cells that produce them, but without bearing the cost of their production. Therefore, it is so important to include in the matrix both impacts produced by the cells itself and by the external factors. Thus, we have an image of the effect of combination therapy on the evolution of cancer cells, which can be a cell killing with recruitment from quiescent (G0) phase or chemotherapy with antiangiotherapy ^[12,31].

In this model the baseline fitness is set to 1, and other parameters used to define the measure of fitness in Tables 6 and 7 are given by:

e cost of producing the cytotoxin

f disadvantage of being affected by the cytotoxin

g benefit of harming other cells

h cost of resistance to the cytotoxin

i cost of proangiogenic factor production

j beneficial effect of receiving the growth factor

r external resources stimulating growth (e.g proangiogenic growth factors)

c external cytotoxic resources (e.g. cytotoxic drugs).

Once more, introduction of resource $r$ leads only to change of the baseline cost of meeting, but resource $c$ effects differently in various measures of fitness. Table 7 contains both resources, $r$ and $c$. The effect of taking into the consideration of proangiogenic external resource will have an effect of changing the baseline fitness from $1$ to $1+r/2$. The alteration of the baseline value normally does not influence the dynamics of the system. However, in this case $r$ parameter can change in time that is why now the baseline will alter in time and space. In terms of tumor cells, the sum of $A$-type and $P$-type may be considered, since $Q$-type does not make any damage to other cells and $R$-type is neutral. On the other hand phenotype A could be considered as cells responsible for the immune system, so then $P$ and $Q$-type shall be tumor cells. In general, the model represents implications of interactions among different cells' phenotypes and feasible stable coexistence. In the Tomlinson model of the so-called angiogenic games, the only condition to reach evolutionary stable state is that the cost of producing growth factors $i$ should be smaller than the benefit $j$. In the extended model (without resources) the expected pay-offs (the sum of the products of frequency and pay-off) are defined as:

To achieve quadruple equilibrium following relations should be satisfied:

Therefore, for a polymorphism (coexistence) of all strategies, each frequency should be contained within the interval (0, 1). It has to be added that calculated formulas for frequencies could be applied only when the above-mentioned conditions are satisfied. In other cases the results could lead to an equilibrium point which may be either an attractor or a repeller, to any other than quadrupled stable polymorphism, to monomorphism or even to unstable populations. To track the evolution of different phenotypes in the population, it is feasible to simulate equations for replicator dynamics ^[16]. They show how frequencies of different strategies change over time, thereby influencing the composition of the studied population. Some examples of the phase portraits (since $A + P + R + Q = 1$, then the graphical representation could be shown as a simplex) are presented on Figures 10 and 11.

For inference analysis in this game, the result when all phenotypes coexist is taken as a reference (Figure 10). Relatively to this result the cost of cytotoxic production has been increased by 0.1 and equals to the benefits from harming the neighbors. Similarly, the adaptation of $P$-cells has been decreased, and at the same time, one of the polymorphic conditions has not been fulfilled. Because of that, $P$-phenotype almost disappears from the population, but the same effect is observed for $Q$-cells. It could be explained by the self-correlation between these two phenotypes – in fact, and the main assumption of the model is that $Q$-cells arise as the evolutionary reaction to the toxic substance produced by $P$-phenotype. It is also observable from the expected results for the quadromorphic population. The fraction of phenotype $P$ is directly proportional to the cost of resistance $h$ and inversely proportional to the losses of interaction with toxic substances $f$. Namely, the more the cells are wounded (including the $P$-cells contact with another $P$-cell and excluding the contact with the resistant $Q$-cells) by the cytotoxic substance, the adjustment of the phenotype $P$ decreases. Similarly with phenotype $Q$, a fraction of which depends on the ratio of the parameters related strictly to phenotype $P$. So the greater the benefits $g$ from the harming of the neighboring cells, the greater the adaptation of $Q$ phenotype. This can be explained by that in contact between $P$ and $Q$, the former does not receive the benefits. On the other hand when parameter $g$ is relatively high (e.g. 0.7), so then $P$ is the dominant strategy in the population displacing all remaining phenotypes (it is the same for the spatial games, which show that when a strategy has got high adaptation, then independently of the type of the game it is a dominant one). Coming back to the results, when $e$ equals $g$ then $P$ and $Q$ cells are displaced from the population, and the game is played between $A$ and $R$ phenotypes in similar terms as in the angiogenic game. Within the reference results, the neutral phenotype $R$ is the dominative one. The analysis could be supported through generating the final frequency of occurrences for parameter changes. This kind of representation does not allow to study the dynamic of phenotype changes in time. However, it is feasible to check the impact of two different parameters for one phenotype at the same time. For example in Figure 12 (the rest of the parameters as in the reference result) an interesting case has arisen for $h = 0$ - phenotype $Q$ decreases while parameter $e$ increases. For $P$-phenotype it can be traced, that if $e$ is greater than 0.4, then $P$ is always 0. It could be related to the fact that value of $g$ in that case equals to 0.4. Increasing $g$ up to 0.7 also gives surprising results.

An important finding is that the four-phenotype model implies third-order dynamics of replication which enables the existence of complex dynamical behaviours including strange attractors. This may be an important hallmark of the evolutionary game theory analysis. A similar model could be used to describe one of the important processes related to communication within cell population resulting from therapy. One very example of such processes analysed by us is the so-called radiation-induced bystander effect. Despite complex analysis of this model due to numerous parameters and relations between cells, the model gives a finite number of diverse results. One of them is the possibility of stable coexistence between different tumor cells within the population. In a similar way, other strategies could be added. However, the graphical representation of the results might be a limitation here. This constraint results from a maximum number of strategies equal to four. What is more, any additional strategies (phenotypes), added to the model, would increase internal dependency between parameters, strategies and results. The machinery of EGT supported by the replication dynamics enables analysis of the evolution of phenotype structure in time within cell populations; nevertheless, it gives no information about the spatial distribution of these phenotypes in tumors. Such possibilities are created by the methodology of spatial evolutionary games theory (SEGT) which enables a study of players' allocation. The lack of perfect mixing is a crucial difference between non-spatial and spatial models. Figure 13 is a spatial counterpart of phase portrait presented on Figure 10 for different reproduction schemes. Within the reference results the neutral phenotype $R$ is the dominative one, then $P$ and also $Q$ appear for the spatial game. There is no obvious explanation for these results since the phenotype $R$ is not better adjusted neither in contact with the rest of phenotypes nor with itself. In the case of the second game (for increased $e$) the phenotype $R$ also dominates for deterministic and quantitative reproductions (the result similar as for non-spatial game). The difference has occurred in the probabilistic reproduction, where $Q$-type has been displaced from the population. Alternatively spatial games could be presented in a way similar to mean-field models. Those outcomes are more focused on the dynamics of the model through the passing generations than on the spatial structures. In that way Figure 14 shows how frequencies of occurrences of each phenotype are changing in time up to the state shown on the Figure 15. What is more the former also confirms the observations and analysis done for the latter and indicates even more clearly the sensitivity with respect to different reproductions. For another analysis, the cost of resistance has been increased by 0.1 compared to the reference model (Figure 14). In this case also one of the polymorphic conditions is not satisfied, and $R$-type is no longer in the population (however in the spatial game it is still dominant phenotype). $P$-cells have increased their frequency of occurrence twice, while $Q$ and $A$ do not change. It shows how difficult is to perform analysis of possible results only by studying the pay-off table without the game simulation (however the number of different parameter sets could also be vast). In this case, increasing $e$ and decreasing $h$ results in regaining the domination by the phenotype $R$ for deterministic and probabilistic reproductions.

The idea of multilayer space evolutionary games which has been presented above allows to take into account heterogeneity of cancer cells. It leads to the conclusion that cancer cells should be considered as representing various phenotypes at the same time described by the frequency of occurrences. This entails more results in the meaning of different spatial structures; however, the analysis could be a bit more difficult. This is why we propose an additional way of graphical representation which takes the average value of players' phenotypes from entire lattice (similar to time charts from the non-spatial game).

Introduction of external resources may result in qualitative changes in the phenotype evolution. In these examples we assume $j = 0$ and $f = 0$. The idea is to not overparametrize the model. Note, that $j$ represents a beneficial effect of receiving the growth factor and $f$ is a disadvantage of being affected by the cytotoxin produced by some cells, and $r$ and $c$ induce similar effects caused by external sources. Thus, the idea is to consider effects of those two factors by changing only two parameters ($r$ and $c$ in this case). The Figures 16-21 present results of the simulations for a different set of parameters.

5.   Discussion

We have proposed a new approach to modelling the effect of interventions on the cancer population evolution based on the evolutionary games with resources. More precisely we apply the multilayer spatial evolutionary games with additional layers representing the evolution of external resources which may describe both the stimulating and the cell-killing treatments. We have posed essentially a question if the genetic heterogeneity in cancer cells exist ab initio or as a result of interventions. Our results, although still far from being complete, indicate that genetic heterogeneity in cancer cells evolve as a result of interventions. However, our response is only partial. More precisely, we can only state that it depends on the character of the intervention. Presented analysis shows that external interventions, depending on its nature are influencing the heterogeneity. One possibility is to change through the change in the baseline of payoff matrix, it means the change of baseline cost of the meeting. The other is an introduction of new phenotypes into the population. This can lead to the conclusion that treatment can generate new phenotypes of cancer cells. This fact means that the initial heterogeneity, that was observed before intervention, in the result of treatment can change producing for example a new, possibly more malignant, cancer subtype. The multilayer approach enables modelling cancer cells as representing different phenotypes at the same time described by a frequency of occurrences. It corresponds to the biologically observed effect of heterogeneity of cancer cells not only at the population level but also at the level of single cells. Moreover, by comparing simulation results for deterministic and probabilistic reproductions, we can study differences in the case of a direct exposure of the cancer cells to external actions and an indirect exposure caused by bystander effect. Mostly the probabilistic reproduction reflects behaviour of the mean field model, as this type of reproduction replace the bystander effect in population individuals. Meanwhile, the deterministic reproduction reflects direct impact of specific subpopulation and results in different model behaviour that the mean field.

To our knowledge, in literature, such discussion is almost absent. An additional key factor (for game theory applications) has been studied by ^[6], being the impact of the ecosystem or the interactions between tumour cells and their environment. They have already modelled changes in the cancer ecosystem in the context of different anti-cancer therapeutic strategies. They indicate that elimination of as many cancer cells as possible may not be essentially the best strategy and have found that destroying only some fraction of the cancer cells (with a particular phenotype) may be far more efficient. Kaznatchev et al. ^[18] is yet another paper in which anticancer treatment scheduling and its effect on the evolution of cancer heterogeneity is considered. The authors focus on the inter-dependence of growth-factor production and acidification and study the resulting mean-field dynamics using the double goods game with three phenotypic strategies. In our study, we are mainly interested in spatial effects, and the number of phenotypes is not an essential constraint. Nevertheless, we also analyze the mean-field case including four phenotypic strategies. Of course, it is an open question as to how to determine a kind of game that cancer cells play and how to experimentally find parameters for game models.

Acknowledgments

This work was supported by the Polish National Science Centre Grants no. DEC-2016/21/B/ST7/02241 (AS), UMO-2015/19/B/ST7/02984 (MKi) the Institute of Automatic Control, Silesian University of Technology under Grant no. BK-204/RAU1/2017, t.3 (MKr) and the National Centre for Research and Development Grant no. STRATEGMED2/267398/4/NCBR/2015 (MILESTONE - Molecular diagnostics and imaging in individualized therapy for breast, thyroid and prostate cancer) (DB).

Calculations were performed on the Ziemowit computer cluster in the Laboratory of Bioinformatics and Computational Biology, created in the EU Innovative Economy Programme POIG.02.01.00-00-166/08 and expanded in the POIG.02.03.01-00-040/13 project.

Conflict of interest

The authors declare there is no conflict of interest.