Evolution of cooperation in an epithelium

Cooperation is prevalent in nature, not only in the context of social interactions within the animal kingdom but also on the cellular level. In cancer, for example, tumour cells can cooperate by producing growth factors. The evolution of cooperation has traditionally been studied for well-mixed populations under the framework of evolutionary game theory, and more recently for structured populations using evolutionary graph theory (EGT). The population structures arising due to cellular arrangement in tissues, however, are dynamic and thus cannot be accurately represented by either of these frameworks. In this work, we compare the conditions for cooperative success in an epithelium modelled using EGT, to those in a mechanical model of an epithelium—the Voronoi tessellation (VT) model. Crucially, in this latter model, cells are able to move, and birth and death are not spatially coupled. We calculate fixation probabilities in the VT model through simulation and an approximate analytic technique and show that this leads to stronger promotion of cooperation in comparison with the EGT model.


Introduction
Tumour development is an evolutionary process whereby cells undergo a series of genetic changes leading to acquired capabilities that confer some growth advantage. In Hanahan and Weinberg's seminal paper [1], six such capabilities or 'hallmarks of cancer' were identified to be necessary for normal cells to become malignant: self-sufficiency in growth signals, insensitivity to antigrowth signals, evading apoptosis, limitless replicative potential, sustained angiogenesis, and tissue invasion and metastasis. Many of these rely on the production of diffusible growth factors [2], the effects of which are felt not only by the producer cell but by other cells in the neighbourhood. As such production of these growth factors can be considered an example of cellular cooperation [3,4].
Models of the evolution of cooperation for diffusible growth factors [5][6][7] have been developed using the framework of evolutionary game theory for well-mixed populations. These models have also been extended to consider spatial effects by placing cells on a lattice [8][9][10] or a fixed graph [11,12]. For the most part these models use periodic boundary conditions; however, there has been some investigation into edge effects [13]. Further examples of the application of game theory to cancer evolution include the reprogramming of energy metabolism [14 -16], micro-environment dependency [17,18], environmental poisoning [19] and invasion [20]. See [21] for a recent review of evolutionary game theory applied to somatic evolution.
Cell populations are not well mixed but organized into tissues or tumours, thus the recent move to incorporate spatial structure is important. Introducing population structure can have a significant effect on evolutionary dynamics [22], in particular, in promoting cooperation [23]. The established framework for modelling games on structured populations, used in the models mentioned above, is evolutionary graph theory (EGT) [24 -30] in which individual cells are placed on the vertices of a graph and neighbours are joined together by edges. Individuals interact and play games with their neighbours, thus deriving their fitnesses. The population evolves via some update rule which dictates how birth and death occur while maintaining the fixed graph structure. When a cell divides, it is necessary for a neighbouring cell to die in order that one of the offspring can occupy the empty vertex. Two commonly used update rules are the birth-death and death-birth rules which essentially differ in the order in which birth and death events occur.
There are several shortcomings of EGT in application to somatic evolution. Tissue and tumour structures are not fixed but dynamic, due to processes of cell division, extrusion and motility. Furthermore, the necessity of births and deaths occurring next to each other is not only unrealistic, but the choice of update rule is one of the main determinants of evolutionary outcomes [31]. Recent work has introduced a new 'shift update' with the aim of addressing the unsuitability of the traditional update rules for cellular structures. The model works extremely well in one dimension [32], predicting enhanced cooperative success compared to other update rules. However, the extension into two dimensions [33] is not straightforward as the shifting of cells disrupts cluster formation of cooperators. This can be resolved by introducing a repulsive force between cells of different types and choosing energy-minimizing shift paths. If the force is strong enough, the shift dynamics is again an effective promoter of cooperation. However, it relies on this somewhat artificial preferential sorting.
Dynamic graph models of evolutionary games also exist; however, they mostly focus on switching connections between vertices, either at random or to increase fitness [34][35][36][37]. These types of models are relevant in social networks, for example, where agents can choose who they interact with and can break social ties with individuals who do not cooperate [38]. They are not good models, however, for populations of cells which are spatially constrained in two-or three-dimensional structures. Furthermore, they still require birth and death to be coupled.
More relevant is the framework developed in [39,40] which uses a topological tissue model [41] to generate a dynamic graph representing cellular interactions. In this dynamical tissue model, birth and death can be spatially decoupled; however, graph topologies do not necessarily correspond to normal tissues. In particular all-defector populations have abnormal polygon distributions, and cooperators on the boundary of a defector cluster can end up with unrealistically high numbers of neighbours. The introduction of forces in a spatial tissue model could resolve these issues.
In order to elucidate what impact, if any, the dynamic nature of cell populations and spatial decoupling of birth and death has on the evolution of cooperation, we will consider evolutionary games on a mechanical model of an epitheliumthe Voronoi tessellation (VT) model [42,43]. Epithelia are the tissues which form the surfaces in the body, such as skin, and the linings of organs. We choose this particular tissue structure as it can be modelled in two dimensions as a sheet of polygonal cells [44]. Furthermore, epithelial cells are highly proliferative compared to other cell types and the source of 85% of cancers making them of particular interest in models of cancer evolution.
We consider cells interacting via an additive Prisoner's Dilemma game, whereby cooperators pay a cost c in order to produce some benefit b for their neighbours. While other games, such as multi-player public goods games, may provide more realistic cancer models, the additive Prisoner's Dilemma is preferred in this study due to its simplicity as a single parameter, two-player game. Furthermore, it is well studied in the EGT context and thus it is straightforward for us to compare with the VT model. In particular, we calculate the fixation probabilities for single mutant cooperators arising in a population of defectors in both models.
We begin, in §2, by introducing EGT and looking at how it can be applied to the evolution of cooperation on epithelia, considering results for an additive Prisoner's Dilemma game with both birth-death and death-birth update rules. We then, in §3, introduce the VT model of an epithelium, again considering the evolution of cooperation under a Prisoner's Dilemma, but this time with spatially decoupled birth and death. We calculate approximate fixation probabilities as well as looking at simulation results. Finally, in §4, we compare these results with the EGT model, finding that cooperation is significantly more successful in the VT model. By running further simulations, implementing an explicit death-birth update in the VT model and a migration analogue into the EGT model, we identify the decoupling of birth and death to be the primary mechanism for the discrepancy.

The model
EGT provides a framework for modelling the evolution of traits on fixed population structures represented by a static graph G. Individuals, labelled i ¼ 1, 2, . . . , N for a population size N, are represented by the vertices, while the edges correspond to neighbour connections. We, therefore, define the adjacency matrix In the additive Prisoner's Dilemma, the trait or type of an individual i is given by s i [ f0, 1g, with s i ¼ 0 denoting a defector (D) and s i ¼ 1 a cooperator (C ). The state of the population is then given by the N-dimensional vector s. For a population in state s, individual i obtains a pay-off f i (s) from its neighbours which is calculated according to a pay-off matrix, given by where b . c and c . 0. The pay-offs are thus the neighbour number). Fitness is then defined to be where d . 0 is the selection strength parameter and the constant 1 takes into account other contributions to fitness. We can let c ¼ 1 without loss of generality, thus the game is defined by a single parameter. Evolution proceeds via a spatial extension of the Moran process [24,45] whereby, at each time step, an individual dies and another reproduces. The offspring occupies the vacant vertex thus keeping the graph structure constant. There are several potential mechanisms for this, known as update rules. Here, we consider two common rules: birth-death: an individual is chosen to reproduce with probability proportional to fitness; its offspring takes the site of a neighbour selected uniformly at random to die; death -birth: an individual is chosen to die uniformly at random; it is replaced by the offspring of a neighbour chosen with probability proportional to fitness.
For a well-mixed population, represented by a complete graph, these two updates rules are equivalent; however, for an arbitrary population structure, the choice of update rule leads to strikingly different dynamics. In the following, we will consider the dynamics in both cases for graph structures representing an epithelium.

Fixation probabilities
In order to consider game dynamics on an epithelium within the EGT context, we consider two different graph structures. Epithelial cells have six neighbours on average; therefore, a hexagonal lattice (HL) is a simple approximation. A VT, however, gives a more realistic representation of an epithelium [46 -48]. There is some variance in neighbour number, but the mean is still 6. The Delaunay triangulation (DT) corresponding to a VT gives the appropriate graph connecting neighbouring cells. See §3.1 and figure 2 for more detail on these terms.
We measure the success of a cooperative mutant by comparing its fixation probability (r C ) to that of a neutral mutant (r 0 ¼ 1/N). Thus if r C . 1/N we say that cooperation is a beneficial mutation or that it is 'favoured by selection'. The critical benefit-to-cost ratio, denoted (b/c)*, is the point where the cooperator fixation probability is equal to the neutral fixation probability, i.e. r C ¼ 1/N.
For a death -birth update rule, we calculate the fixation probabilities against benefit-to-cost ratio (b/c) for an HL and DT with a population size of N ¼ 100 and periodic boundary conditions. Results are plotted in figure 1 in which each data point is the result of 1 Â 10 5 simulations.
Analytical results are calculated using the theory developed in [30], where the authors derive an equation for the fixation probabilities on any graph. Here, t n is the expected coalescence time from the two ends of an n-step random walk, where the initial vertex is chosen proportional to degree. Thus these quantities are purely properties of the graph and can be calculated computationally by solving a recurrence relation. We use a small selection strength, d ¼ 0.025, and there is a good fit between simulation and theory in the range shown for b . 4. Furthermore, the heterogeneity in the DT seems to have a negligible effect on fixation probabilities compared to the dependence on benefit-to-cost ratio.
Critical benefit-to-cost ratios are calculated for both graphs from simulations and equation (2.5) and summarized in table 1.
The results are very different for a birth-death update rule: cooperation is never favoured by selection under an additive Prisoner's Dilemma game and r C , 1/N for all b , c, c . 0 [25,28,31]. Thus within the EGT framework cooperation is only a successful evolutionary strategy on an epithelial structure with a death-birth update above a critical benefit-to-cost ratio of approximately 6.7.
The HL seems to be a reasonable approximation to the structure. Using the more realistic DT with neighbour number  heterogeneity does not significantly alter fixation probabilities or the critical benefit-to-cost ratio, at least in the weak selection limit we are using. We note, however, that these results are for an average pay-off and that an accumulative pay-off (in which pay-offs are simply summed over interactions) can amplify differences due to heterogeneity. We should also note that cooperation is possible in well-mixed populations or graphstructured populations with birth-death update for games other than the Prisoner's Dilemma, such as the snowdrift or stag-hunt games, and it is possible to generalize (2.5) to analyse these [30]. Whether or not these results are illuminating in terms of a real epithelium is an important question, however, and as we have noted previously there are some serious shortcomings to the model, first that population structure is static and secondly the troubling dependence on the update rule. Which update rule is closest to reality is unclear and while there likely is some coupling in birth and death processes in a real epithelium, there is certainly no absolute requirement for birth and death events to occur next to each other. In order to explore whether these factors are important to the dynamics we will move on to consider the VT model of an epithelium in which cells are able to move past each other and birth and death are spatially decoupled.

Voronoi tessellation model of an epithelium
In order to analyse the dynamics of evolutionary games on a more realistic population structure, we will use the VT model [42,43] developed for the colonic crypt epithelium. In the following, we will explain how the mechanical model works and generates a time-dependent graph structure on which to study evolutionary game dynamics. We will then derive an approximation for the fixation probability and use these results along with simulation to compare with the EGT model.

The model
The VT model represents a tissue as a set of points corresponding to the centres of individual cells. These points lie in a fixed domain with periodic boundary conditions. Cells move freely in space and exert spring-like forces on one another, such that is the force exerted by cell j on its neighbour i. Here, m is the spring constant and r ij ¼ r i 2 r j , where r i is the position vector of cell i andr ij is the corresponding unit vector. The natural separation between cells s ij (t) ¼ s is constant and the same for all neighbour pairs. The exception to this is for newborn sister cells for whom s ij grows linearly from e to s over the course of an hour. The total force acting on cell i is then where N i (t) is the set of cells neighbouring i. By assuming that motion is over-damped due to high levels of friction we obtain the equation of motion for each cell in the form of a first-order differential equation where h is the damping constant. This is solved numerically using where Dt is a sufficiently small time step for numerical stability. Parameter values used in our simulations are taken from [44] and based on studies of the colonic crypt [42,49]. These are summarized in table 2. While changes in these values affect the dynamics, our main result is robust (see the electronic supplementary material information where we consider changes to m). The neighbour connections between cells are determined by the VT of the set of cell-centres (figure 2). The VT divides the plane into polygons, where each polygon is defined as the region of the plane closer to its generator (i.e. cell-centre) than any other. Each cell can, therefore, be represented as a distinct region with a well-defined area and neighbour set. The dual graph to the VT is the DT in which the cell centres are the graph vertices and neighbours are connected by edges. The DT, therefore, gives the adjacency matrix A ij (t) from which we can calculate cell fitnesses. As it is defined by the cell-centre positions, the DT must be recalculated after every time step during which cells may have moved, died or reproduced.
As in the previous model, we allow the system to evolve by a Moran process whereby birth and death events occur simultaneously. The key difference is that we decouple the locations of these events. We also implement the process in continuous rather than discrete time, noting that a translation to continuous time in the previous model does not affect fixation probabilities [30] and therefore the results are directly comparable. In the continuous time Moran process, update events occur at exponentially distributed times with rate l. When an update event occurs, a mother cell is chosen at random from the population with probability proportional to fitness. This cell divides creating two daughter cells, which are exact clones of the mother. A cell is also chosen to die (i.e. to be extruded from the tissue) uniformly at random. This process is represented in figure 3.
To calculate fixation probabilities for a single mutant cooperator invading a defector population in the VT model, we run simulations as follows. We begin with defector cells placed on a regular HL with periodic boundary conditions and the simulation algorithm proceeds until the system has relaxed into a dynamic equilibrium. We then choose a random cell to become a cooperator and continue the simulation until only cooperators or defectors remain. The simulation algorithm consists of the following steps: (i) DT is performed to determine

Approximating the fixation probabilities
Due to the complexities of the VT model, it is not possible to derive exact analytical solutions as was done for EGT [30]. Instead, we look for approximate solutions by considering the expected fitness for different cell types in populations with a given number of cooperators [50]. While the graph is dynamic and dependent on the spatial distribution of points, it is also planar and mechanically constrained. Furthermore, if we begin with a single mutated cell, its progeny are likely to remain in a cluster as the clone grows. Thus we assume that variation in fitnesses for cells of each type will be small for a given number of cooperators in the population and that the average over a large number of states is a good approximation. Comparing our theoretical results to simulations, we find that fixation probabilities calculated based on this assumption are good approximations. Let us denote a state with n cooperators S n ¼ (s n , G), where s n is the vector of cell types and G is the graph. Then we define T þ/2 (S n ) to be the probability that when an event occurs the number of cooperators is increased/decreased by one, i.e. and We can then define the average transition probabilities for a state with n cooperators to be T + n ¼ hT + (S n )i where the average is taken over a large ensemble of possible states. Substituting in for the fitnesses (2.4) and taking the weak selection limit d ( 1 we obtain and where k.l 0 denotes an average over a large ensemble of possible states for the neutral process d ¼ 0 and are the average cooperator fitness and average fitness, respectively. From (2.3) and (3.9), we obtain is the normalized average number of degree-weighted cooperator-cooperator interactions in a system with n cooperators. This can be calculated computationally by running simulations for a neutral process and tracking clones (groups of cells with common ancestry). At each time interval, we calculate the contribution to L CC n for all clones in the system, treating each lineage as a group of n cooperators in a population of defectors. See figure figure 4 for a plot of L CC n with N ¼ 100.
We use the equation for cooperator fixation probability derived in [51] for a well-mixed population (3:12) Figure 3. Spatially decoupled update rule in the Voronoi tessellation model. When an update event occurs, a mother cell is chosen to reproduce with probability proportional to fitness (blue). A second cell is chosen to die uniformly at random (red). The mother cell divides and the dead cell is removed from the tissue. (Online version in colour.)  royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 16: 20180918 In that case, the transition probabilities and thus g n are defined exactly for each value of n. For the VT model, we substitute in the mean transition probabilities given by equations (3.7), (3.8) and (3.9), to obtain for the fixation probability in the weak selection limit. The critical benefit-to-cost ratio is then obtained by setting : (3:14) Figure 5 compares equation (3.13) with simulation results for the VT model. It shows there is a reasonable fit between our approximation of fixation probabilities with the simulation data in the region 2.0 , b , 3.5, where we have once again set c ¼ 1. These values are close to the critical benefit-to-cost ratio and therefore represent the region in which we would expect the weak selection limit to hold, thus this equation for fixation probabilities is a reasonable approximation. The critical benefit-to-cost ratios calculated from simulation and equation (3.14) are given in table 1. For both, we get a value of b/c ¼ 2.8 correct to one decimal place. This is significantly less than the critical benefit-to-cost ratios calculated for the EGT model with death-birth update. In the next section, we will look further at comparing these models and attempt to identify the mechanism by which cooperation is promoted in the VT model. Figure 6 shows the results of these simulations along with the theoretical EGT results for the HL graph with death-birth update and the critical benefit-to-cost ratios are summarized in table 1. It is clear that cooperators are much more successful in the VT model, in particular, the critical benefit-to-cost ratio for the VT model is less than half that for EGT with death-birth update. The question then arises as to what mechanism is causing this amplifying effect in the VT model, the two obvious candidates being the effect of cell motility and the decoupling of birth and death. One way to test whether cell motility is enhancing the evolutionary success of cooperation is to introduce an analogue into the EGT model whereby we allow cells to swap sites with their neighbours. At each time step, a swap occurs with probability m. When this happens, a cell is chosen uniformly at random to switch places with one of its neighbours. Note that this process is independent of cell fitness. The parameter m is, therefore, a measure of the strength of migration and by setting m ¼ 0 we regain the original EGT model. Figure 7 plots fixation probability against  Figure 5. An approximation for fixation probabilities in the VT model is given by equation (3.13) and plotted here (solid line) for d ¼ 0.025, c ¼ 1 and N ¼ 100. Comparison with simulation results ( points) shows that the approximation is good near the critical benefit-to-cost ratio (i.e. where r C ¼ r 0 ¼ 1/N), but breaks down outside the region 2 , b , 3.5. This is consistent with the fact that the equation is derived in the weak selection limit, and suggests that it can be used to calculate the critical ratio. (Online version in colour.)  Figure 6. Fixation probabilities for the Prisoner's Dilemma game in the VT model with c ¼ 1 and d ¼ 0.025. Points show simulation results for a decoupled update rule (blue, circles) and a death-birth update rule (red, squares). For the decoupled update rule, the approximate fixation probabilities given by equation (3.13) are plotted (blue, solid line) and for the death-birth update we plot a best fit line (red, dashed line). Fixation probabilities, given by equation (2.5), for an HL with death-birth update in the EGT model (green, solid line) are also shown for comparison. The grey dotted line shows the fixation probability for a neutral mutant. It is clear that cooperation is significantly favoured in the VT model with decoupled update rule when compared with the EGT results, in particular, the critical benefit-to-cost ratio is more than halved. However, when a death-birth update is introduced on the VT model this effect is suppressed and the critical benefitto-cost ratio is very close to the EGT case. (Online version in colour.)  benefit-to-cost ratio for a range of m values and demonstrates that increasing migration within this framework actually decreases the evolutionary success of cooperation. It, therefore, seems unlikely that the ability of cells to move past each other in the VT model is the reason for enhanced cooperative success.

Comparing the models
In order to determine whether the spatial decoupling of birth and death promotes cooperation, we consider the VT model with a death-birth update rule. To implement this, we follow the simulation algorithm as defined in §3, the only change being in choosing which cells reproduce and die when an update event occurs. First, a cell is chosen for extrusion uniformly at random. Fitnesses are then calculated for the neighbouring cells and one of these is chosen to divide with probability proportional to fitness. This process is shown schematically in figure 8. It can be seen clearly in figure 6 that changing the update rule in this way suppresses the evolutionary success of cooperation in comparison to the decoupled update rule. Indeed, in this case, we obtain b/c ¼ 7.3 which is greater than for the EGT model with death-birth update.
Combining these two results, we conclude that it is the spatial decoupling of birth and death which leads to the amplification of cooperative success in the VT model. Indeed, this is an intuitive result and is consistent with results from the shift dynamics models [32,33]. A cooperative strategy is only beneficial if cells are able to form a cluster of cooperators. If birth and death are constrained to occur next to each other, as is the case for death-birth and birth-death update rules, then the cluster can only grow at the boundary. If a cell were to reproduce inside a cooperative cluster, it would result in the death of a neighbouring cooperator, leaving the size of the cooperator population unchanged. For the decoupled birth and death update in the VT model, this is not the case. If a cooperator inside the cluster reproduces, it will lead to an increase in the size of the cooperator population with probability 1 2 n/N, where n is the number of cooperators and N the total number of cells. The fact that migration appears to suppress the success of cooperation could also provide an explanation as to why, if a death-birth update is enforced in both cases, cooperators fare better in the EGT model than in the VT model.

Conclusion
EGT has become the accepted framework for modelling the evolution of cooperation on structured populations, ranging from complex social networks to collective cellular behaviour organized in tissues. While it may be an appropriate tool for the former, we have demonstrated that a static graph model is not sufficient to capture the dynamic behaviour of an epithelium.
We have shown using the theory developed by Allen et al. [30] and simulations that for a Prisoner's Dilemma on an epithelium-like structure in EGT, cooperation is successful if b/c . 6.7 for a death-birth update, where we have used an averaged pay-off. This inequality holds when we model the epithelium as an HL as well as a DT, suggesting that there is a marginal effect on fixation probabilities due to heterogeneity of neighbour number. However, the choice of an averaged payoff could be suppressing the effect of heterogeneity compared to an accumulated pay-off, as it does for scale-free networks [52,53]. It would be advisable therefore to compare fixation probabilities on the two structures for an accumulated payoff, although we do not expect a substantial difference. Vertex degree in scale-free networks follow a polynomial distribution and therefore exhibit large variance, whereas degree variance in DTs is comparatively small.
For a birth-death update on the other hand, cooperation is not successful for any benefit-to-cost ratio under a Prisoner's Dilemma game. The fact that the dynamics are so sensitive to the choice of update rule is troubling and neither update rule is a realistic representation of birth and death in an epithelium. For the VT model, we are able to spatially decouple birth and death. We showed, using simulation and approximate theoretical results, that using a decoupled update rule in the VT model promotes cooperation compared to the EGT examples. Furthermore, when the VT model was run with a death-birth update this effect was suppressed and cooperation actually fared worse than in the EGT model, leading us to conclude that the decoupling of birth and death is the main mechanism for increased success of cooperation in the VT model. This is consistent with previous work looking at shift dynamics on a static graph which found that decoupling birth and death led to increased cooperative success in one dimension [32], and in two dimensions if a repulsive force was introduced between cells of different types [33]. The fact that cells can move and change neighbours in the VT model, however, does not appear to increase the likelihood of cooperation fixating. Indeed we found that introducing migration into an EGT model actually suppressed cooperation, and it is, therefore, possible that cell motility is acting to reduce cooperative success in the VT model.
As it is the update rule which seems to influence the evolutionary success of cooperation most substantively, the question arises as to which, if any, reflects the behaviour of a real epithelium. Clearly, it is unrealistic that when a death Figure 8. Death -birth update rule in the Voronoi tessellation model. When an update event occurs, a cell is chosen to die uniformly at random from the population (red). From the neighbourhood of the dead cell (yellow), a mother cell (blue) is then chosen with probability proportional to fitness. The mother cell divides and the dead cell is removed from the tissue. (Online version in colour.) royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 16: 20180918 occurs, it is immediately followed by a neighbour undergoing division, or vice versa, as for the death-birth and birth-death update rules respectively. However, it is also not the case that birth and death events are completely spatially independent. Cell extrusion can be induced in areas of overcrowding within a tissue, which could be caused by high levels of proliferation. Similarly, if local density is low, e.g. due to a high instance of cell death, cells can be induced to reproduce [54,55]. It is difficult to see how this more subtle link between birth and death could be implemented in an EGT model; however, the VT model could be extended to include density-dependence for division and/or extrusion. Furthermore, a density-dependent model would allow us to maintain an (almost) constant population size without enforcing that birth and death occur simultaneously, another unrealistic assumption.
In our discussion of whether cooperation is successful on an epithelium we have limited ourselves to the additive Prisoner's Dilemma game, whereas evolutionary game theory models of cancer have used a variety of social dilemma games. Extending our analysis to a general two-strategy game should be relatively straightforward, indeed we can use the critical benefit-to-cost ratio to calculate the structure coefficient and derive a general condition for evolutionary success for a two-player, two-strategy game [56]. However, it has been argued that multi-player public goods games are more realistic for cancer modelling, and can lead to very different results. Recent work has considered the dynamics of these types of games on lattices [9] and DT graphs [12] in an EGT framework, it would, therefore, be an interesting comparison, but non-trivial extension, to consider them on the VT model. Data accessibility. The code and data can be accessed at https://github. com/jessiesrr/evo-epithelium.
Authors' contributions. J.R. and K.M.P. designed the research. J.R. carried out the research and wrote the paper. K.M.P. edited the paper.
Competing interests. We declare we have no competing interests. Funding. This research was funded by an EPSRC studentship held by J.R.