Evolution of cooperation on an epithelium

TL;DR

This paper compares the evolution of cooperation in epithelial tissue models, showing that cell mobility enhances cooperative success more than static population models, through simulations and analytical methods.

Contribution

It introduces a dynamic epithelial model using Voronoi tessellation and demonstrates that cell movement promotes cooperation more effectively than traditional static models.

Findings

Cell mobility in tissue models enhances cooperation.

Voronoi tessellation model predicts higher fixation probabilities.

Dynamic models outperform static graph models in cooperation promotion.

Abstract

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. 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 evolutionary graph theory, to those in a mechanical model of an epithelium =- the Voronoi tessellation 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 Voronoi tessellation model through simulation and an approximate analytic technique and show that this leads to stronger promotion of cooperation in comparison with the evolutionary graph theory model.

Figures

Tables

Table 1. Table 1 : Summary of critical benefit-to-cost ratios, ( b / c ) ∗ superscript 𝑏 𝑐 (b/c)^{*} , for the different models: a Moran process with death-birth update on a Voronoi network and a hexagonal lattice; a Voronoi tessellation model. Results are shown for both the theory and simulations.

	Theory	Simulation
EGT model with DT (death-birth)	6.69	6.74
EGT model with HL (death-birth)	6.68	6.67
VT model (decoupled update)	2.78	2.83
VT model (death-birth)	-	7.26

Table 2. Table 2 : Table of parameters used in the Voronoi tessellation model [ 42 ] .

Parameter	Description	Value
$\mu$	Spring constant	50.0
$s$	Natural seperation of mature cells	1.0
$ϵ$	Initial seperation of sister cells	0.05
$\eta$	Drag coefficient	1.0
$\Delta t$	Time step (h)	0.005
$\lambda$	Division and apoptosis rate (h-1)	${12.0}^{-1}$

Equations

$$A_{ij}=\{\begin{array}{ll}1,& \text{ if i and j are neighbours}\\ 0,& \text{ otherwise.}\end{array}$$

$$\blockarray{ccc}&CD\\ \block{c(cc)}Cb-c-c\\ Db0\\ \;,$$

$$f_i(\mathbf{s})=-c{s}_i+b\sum _{j\in G}\frac{A_{ij}{s}_j}{k_i},$$

$$F_i(\mathbf{s})=1+\delta f_i(\mathbf{s}),$$

$${\rho }_C=\frac{1}{N}+\frac{\delta }{2N}\left(-c{t}_2+b(t_3-t_1)\right)+\mathcal{O}({\delta }^2)$$

$${\mathbf{F}}_{ij}(t)=\mu {\hat{\mathbf{r}}}_{ij}(t)(|{\mathbf{r}}_{ij}(t)|-s_{ij}(t))$$

$${\mathbf{F}}_i(t)=\sum _{j\in {\mathcal{N}}_i(t)}{\mathbf{F}}_{ij},$$

$$\eta \frac{d{\mathbf{r}}_i}{dt}={\mathbf{F}}_i(t),$$

$${\mathbf{r}}_i(t+\Delta t)={\mathbf{r}}_i(t)+\frac{\Delta t}{\eta }{\mathbf{F}}_i,$$

$$T^{+}(S_n)$$

$$T^{-}(S_n)$$

$$T_n^{+}$$

$$T_n^{-}$$

$$f_C=\frac{1}{n}\sum _{i\in G}{s}_i{f}_i$$

$$⟨f_C-f{⟩}_0=-c\left(1-\frac{n}{N}\right)+b\left({\Lambda }_n^{CC}-\frac{n}{N}\right),$$

$${\Lambda }_n^{CC}=\frac{1}{n}{⟨\sum _{i,j\in G}\frac{s_i{s}_j{A}_{ij}}{k_i}⟩}_0$$

$${\rho }_C=\left[1+\sum _{m=1}^{N-1}\prod _{n=1}^m{\gamma }_n\right]$$

$${\rho }_C\approx \frac{1}{N}+\frac{\delta }{N}\left\{\frac{-c(N-1)}{2}+b\sum _{m=1}^{N-1}\sum _{n=1}^m\left(\frac{{\Lambda }_n^{CC}-n/N}{N-n}\right)\right\}+\mathcal{O}({\delta }^2)$$

$${\left(\frac{b}{c}\right)}^{\ast }\approx \frac{N-1}{2}{\left[\sum _{m=1}^{N-1}\sum _{n=1}^m\left(\frac{{\Lambda }_n^{CC}-n/N}{N-n}\right)\right]}^{-1}.$$

Full text

**Keywords— ** cooperation, epithelium, population structure, evolutionary graph theory, Voronoi tessellation

1 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].

The reprogramming of energy metabolism is also considered to be a hallmark of cancer [5] and refers to the fact that cancer cells tend to metabolise through glycolysis rather than aerobic respiration, even when oxygen is abundant. This is known as the Warburg effect [6]. It has been postulated that glycolytic cells can be considered as cooperators, in that they produce lactic acid as a by-product which is toxic to healthy cells, and can thus be considered a shared benefit to the cancer cells.

Models of the evolution of cooperation for diffusible growth factors [7, 8, 9] and the Warburg effect [10, 11, 12] 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 [Bach, 14, 15] or a fixed graph [16, 17]. Further examples of the application of game theory to cancer evolution include microenvironment dependency [18, 19], environmental poisoning [20] and invasion [21]. See [22] for a recent review of evolutionary game theory applied to somatic evolution.

Cell populations are not well-mixed but organised 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 [23], in particular in promoting cooperation [24]. The established framework for modelling games on structured populations, used in the models mentioned above, is evolutionary graph theory (EGT) [25, 26, 27, 28, 29, 30, 31] 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 evolutionary graph theory 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 [32]. 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 [33], predicting enhanced cooperative success compared to other update rules. However the extension into two-dimensions [34] 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-minimising 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 focus on switching connections between vertices, either at random or to increase fitness [35, 36, 37, 38]. 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 [39]. 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.

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 epithelium- the Voronoi tessellation (VT) model [40, 41]. 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 [42]. 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.

Rather than focussing on a particular cancer model, we consider the simple, and well studied, example of an additive prisoner’s dilemma game, whereby cooperators pay a cost $c$ in order to produce some benefit $b$ for their neighbours. We consider whether the results for the VT model are significantly different from those obtained from EGT. In particular we calculate the fixation probabilities for single mutant cooperators arising in a population of defectors in both models.

We begin, in Section 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 Section 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 Section 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.

2 Evolutionary graph theory

2.1 The model

Evolutionary graph theory 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

[TABLE]

In the additive prisoner’s dilemma, the trait or type of an individual $i$ is given by $s_{i}\in\{0,1\}$, 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 $\mathbf{s}$.

For a population in state $\mathbf{s}$, individual $i$ obtains a payoff $f_{i}(\mathbf{s})$ from its neighbours which is calculated according to a payoff matrix, given by

[TABLE]

where $b>c$ and $c>0$. The payoffs are thus

[TABLE]

where $k_{i}=\sum_{j\in G}A_{ij}$ is the degree of vertex $i$ (i.e. the neighbour number). Fitness is then defined to be

[TABLE]

where $\delta>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 [43, 25] 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.

2.2 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 Voronoi tessellation however gives a more realistic representation of an epithelium [44, 45, 46]. 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 Section 3.1 and Figure 2 for more detail on these terms.

We measure the success of a cooperative mutant by comparing its fixation probability ($\rho_{C}$) to that of a neutral mutant ($\rho_{0}=1/N$). Thus if $\rho_{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. $\rho_{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\text{\times}{10}^{5}$ simulations. Analytical results are calculated using the theory developed in [31], where the authors derive an equation

[TABLE]

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, $\delta=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. These are calculated for both graphs from simulations and Equation (5) and summarised 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 $\rho_{C}<1/N$ for all $b<c$, $c>0$ [26, 29, 32]. 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 payoff and that an accumulative payoff (in which payoffs are simply summed over interactions) can amplify differences due to heterogeneity. We should also note that cooperation is possible in well-mixed populations or graph structured 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 generalise (5) to analyse these [31].

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.

3 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 [40, 41] 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.

3.1 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

[TABLE]

is the force exerted by cell $j$ on its neighbour $i$. Here $\mu$ is the spring constant and $\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$, where $\mathbf{r}_{i}$ is the position vector of cell $i$ and $\hat{\mathbf{r}}_{ij}$ is the corresponding unit vector. The natural seperation 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 $\epsilon$ to $s$ over the course of an hour.

The total force acting on cell $i$ is then

[TABLE]

where $\mathcal{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

[TABLE]

where $\eta$ is the damping constant. This is solved numerically using

[TABLE]

where $\Delta t$ is a sufficiently small time step for numerical stability. For the parameter values used in our simulations see Table 2.

The neighbour connections between cells are determined by the VT of the set of cell-centres (see Figure 2). The VT divides the plane into polygons, where each polygon is defined as the 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 Delaunay triangulation (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 timestep 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 [31] and therefore the results are directly comparable. In the continuous time Moran process update events occur at exponentially distributed times with rate $\lambda$. 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 hexagonal lattice 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: (1) DT is performed to determine cell neighbours; (2) forces are calculated and the cells moved accordingly; (3) an update event occurs with probability $N\lambda\Delta t$, in which case fitnesses are calculated according to the evolutionary game and the decoupled update rule is applied.

3.2 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 [31]. Instead we look for approximate solutions by considering the expected fitness for different cell types in populations with a given number of cooperators [47]. 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}=(\mathbf{s}_{n},G)$, where $\mathbf{s}_{n}$ is the vector of cell types and $G$ is the graph. Then we define $T^{+/-}(S_{n})$ to be the probability that when an event occurs the number of cooperators is increased/decreased by one, i.e

[TABLE]

We can then define the average transition probabilities for a state with $n$ cooperators to be $T^{\pm}_{n}=\langle T^{\pm}(S_{n})\rangle$ where the average is taken over a large ensemble of possible states. Substituting in for the fitnesses (4) and taking the weak selection limit $\delta\ll 1$ we obtain

[TABLE]

where $\langle.\rangle_{0}$ denotes an average over a large ensemble of possible states for the neutral process $\delta=0$ and

[TABLE]

are the average cooperator fitness and average fitness respectively. From (3) and (14) we obtain

[TABLE]

where

[TABLE]

is the normalised 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 $\Lambda^{CC}_{n}$ for all clones in the system, treating each lineage as a group of $n$ cooperators in a population of defectors. See Figure 4 for a plot of $\Lambda^{CC}_{n}$ with $N=100$.

We use the equation for cooperator fixation probability derived in [48] for a well-mixed population

[TABLE]

with $\gamma_{n}=T^{-}_{n}/T^{+}_{n}$. In that case the transition probabilities and thus $\gamma_{n}$ are defined exactly for each value of $n$. For the VT model we substitute in the mean transition probabilities given by Equations (12), (13) and (14), to obtain

[TABLE]

for the fixation probability in the weak selection limit. The critical benefit-to-cost ratio is then obtained by setting $\rho_{C}=1/N$ giving

[TABLE]

Figure 5 compares Equation (18) 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 (19) 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.

4 Comparing the models

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 summarised 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 timestep 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 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.

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 Section 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 [33, 34]. 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-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.

5 Conclusions

Evolutionary graph theory has become the accepted framework for modelling the evolution of cooperation on structured populations, ranging from complex social networks to collective cellular behaviour organised 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. [31] 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 payoff. 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 payoff, as it does for scale-free networks [49, 50]. 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 is so sensitive to the choice of update rule is troubling and neither update rule is a realistic representation of birth and death on 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 [33], and in two-dimensions if a repulsive force was introduced between cells of different types [34]. 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 cooporation, 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 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 [51, 52]. 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 succesful 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 [53]. However it has been argued that multiplayer 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 [14] and DT graphs [17] 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. JR and KMP designed the research. JR carried out the research and wrote the paper. KMP edited the paper.

Acknowledgements. We would like to thank John Talbot for helpful discussions and Pilar Guerrero and David Page for advice on the code.

Funding. This research was funded by an EPSRC studentship held by JR.

Competing interests. We declare we have no competing interests.